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AFIT/GE/EN  G /92D-37 


Abstract 


The  Air  Force  has  equipped  its  aircraft  with  avionic  systems  such  as  Globa, 1  Posi¬ 
tioning  Systems  (GPS)  and  Inertial  Guidance  Systems  (INS)  capable  of  providing  accurate 
navigation  solutions.  The  aircrews  hying  these  aircraft  require  a  system  that  can  either 
survive  the  hostile  environments  encountered  in  combat  or  notify  the  aircrew  that  their 
performance  has  been  significantly  degraded.  This  research  focuses  on  failure  detection  and 
isolation  techniques  using  u.n  extended  Kalman  filter  and  generalized  likelihood  ratios  using 
matched  filters.  Analysis  is  conducted  using  a  Kalman  filter  development  package  known 
as  the  Multimode  Simulation  for  Optimal  Filter  Evaluation  (MSOFE).  Both  a  large  order 
truth  model  for  the  navigation  system  (in  which  a  full  24  satellite  constellations  is  modeled) 
and  a  reduced-order  Kalman  filter  are  developed.  Results  suggest  that  failures  within  the 
GPS  can  be  detected,  isolated,  and  in  some  cases  compensated  through  feedback. 
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DETECTION  OF  SPOOFING,  JAMMING,  OR  FAILURE 
OF  A  GLOBAL  POSITIONING  SYSTEM  (GPS) 


1.  Introduction 


A  variety  of  Global  Positioning  System  (GPS)  receivers  and  Inertial  Navigation  Sys¬ 
tems  (INS)  are  installed  on  military  aircraft.  The  GPS  receives  information  from  orbiting 
satellites  and  calculates  estimates  of  the  position  and  velocity  of  the  aircraft.  The  INS 
detects  inertial  motion  of  the  aircraft  and  calculates  its  own  estimates  of  aircraft  position 
and  velocity.  The  GPS  and  INS  are  integrated  to  form  a  system  that  is  more  accurate  and 
reliable  than  each  of  these  systems  by  itself.  Additional  measurements  from  ground  based 
uiauspoiiueis  are  avanaule  ami  am  in  de^erming  aircraft  position.  The  transponder  based 
system  is  referred  to  as  the  Range/ Range  -Rate  System  (RRS).  The  combination  of  the 
GPS,  INS,  and  RRS  form  the  Navigation  Reference  System  (NRS)  whose  primary  function 
is  to  assist  the  pilot  in  navigating  the  aircraft.  An  important  note  is  that  the  RRS  is  often 
used  on  test  ranges  to  reconstruct  flight  paths  over  the  range.  It  is  not  desired  to  limit  the 
applications  of  this  thesis  to  test  aircraft,  so  it  is  assumed  that  the  RRS  represents  any 
one  of  many  ground-based  transponder  systems  used  on  current  military  aircraft. 

1.1  Background 


Research  is  being  done  by  the  Central  Inertial  Guidance  Test  Facility  (CIGTF), 
6585th  Test  Group,  Air  Force  Systems  Command  (AFSC),  Holloman  AFB,  NM  to  deter¬ 
mine  the  vulnerability  of  the  GPS  to  jamming  and  spoofing.  Jamming  is  nothing  more 
than  bombarding  the  GPS  receiver  with  electronic  noise.  Figure  1.1  illustrates  jamming 
as  a  broad  band  of  noise  directed  at  the  GPS  receiver.  There  is  no  intent  to  mislead  the 
GPS  receiver  with  jamming,  but  rather  to  prevent  it  from  receiving  the  desired  satellite 
information.  Spoofing  is  more  complex  and  harder  to  detect.  The  goal  of  spoofing  is  to 
mimic  the  signal  sent  from  the  satellite  to  the  GPS  receiver  but  with  minor  changes  to  the 
signal.  These  minor  changes  will  cause  the  GPS  to  calculate  erroneous  estimates  and  draw 
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the  aircraft  away  from  its  desired  destination.  Spoofing  is  shown  in  Figure  1 .1  as  a  directed 
signal  from  a  ground -based  platform  aimed  at  a  particular  target.  An  airborne  platform 
for  the  spoofer  would  be  more  effective  against  aircraft  with  directional  GPS  antennas 
capable  of  pointing  away  from  a  ground-based  spoofer.  Jamming  and  spoofing  constitute 
two  of  the  three  failures  that  will  be  addressed  in  this  thesis.  A  third  failure  that  can  affect 
the  performance  of  the  GPS  is  the  lose  of  a  pseudorange  input  to  the  receiver.  Failures 
caused  by  losing  a  GPS  pseudorange  may  look  similar  to  spoofing  or  jamming,  but  it  is 
important  to  identify  which  of  the  three  failures  has  occurred. 


Enemy 


1.2  Pwblem  Definition 

The  primary  goal  of  this  thesis  will  be  to  develop  a  failure  detection  and  isolation 
(FDI)  system  that  will  identify  GPS  failures  as  described  above.  This  FDI  system  will 
be  based  on  techniques  discussed  in  the  literature  review.  The  follow-on  goal  will  be  to 
develop  an  adaptive  system  which  can  correct  for  a  failure  that  has  been  delected  and 
isolated.  The  NR.S  is  to  retain  good  performance  as  a  result  of  corrections  which  are 
typically  used  in  a  feedback  configuration. 
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1.3  Scope 


Two  types  of  failures  will  be  considered.  Figure  1.2  depicts  the  two  types  of  fail¬ 
ures  as  a  jump  or  a  ramp.  A  jump  failure  results  in  a  corrupted  signal  with  a  different 
magnitude  than  the  desired  signal  and  occurs  almost  instantaneously  in  time.  Some  FDI 
algorithms  assume  the  magnitude  of  this  jump  is  known  but  it  is  more  realistic  to  assume 
the  magnitude  is  unknown  and  important  for  the  FDI  system  to  determine.  A  ramp  fail¬ 
ure  increases  more  gradually  than  a  jump.  The  deceiving  signal  is  constantly  increasing 
away  from  the  desired  signal  at  an  unknown  rate.  The  thresholds  and  delay  shown  will  be 
discussed  in  the  literature  review.  The  rate  of  increase  of  the  ramp  and  magnitude  of  the 
jump  will  be  varied  for  different  studies. 

Although  the  failure  algorithms  discussed  in  this  thesis  could  apply  to  systems  other 
than  the  GPS,  no  attempt  will  be  made  to  induce  failures  into  the  INS  or  RRS,  as  the 
primary  focus  will  be  on  the  GPS  only.  Criteria  for  completion  of  the  research  are  presented 
in  Section  1.6. 
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Figure  1.2.  Types  of  Failures 
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1-4  Assumptions  and  System  Description 

The  list  that  follows  helps  describe  the  scenario  chosen  for  this  thesis  and  the  simpli¬ 
fying  assumptions  that  were  adopted.  The  impact  caused  by  these  assumptions  may  not 
be  apparent  without  knowing  the  context  in  which  they  were  applied,  so  references  arc 
listed  to  help  the  reader  see  their  significance. 

1.  The  NRS  will  be  mounted  on  a  computer-simulated  aircraft  that  is  capable  of  high 
dynamic  maneuvers  analogous  to  a  military  fighter  jet.  Actual  flight  tests  will  not 
be  performed,  but  a  two-liour  flight  profile  as  shown  in  Figure  C.l  will  be  simulated 
on  a  computer.  This  profile  is  represented  by  the  latitude,  longitude,  and  altitude  of 
the  aircraft. 

2.  A  typical  NRS  configuration  will  be  used  as  shown  by  a  block  diagram  representa¬ 
tion  in  Figure  1.3.  Error  states  are  used  as  the  basis  of  the  filter  design  model,  and 
difference  measurements  are  provided  to  a  Kalman  filter  by  the  GPS.  INS,  and  RRS. 
The  Kalman  filter  generates  error  state  estimates  used  to  correct  the  original  INS 
states,  resulting  in  refined  estimates  of  latitude,  longitude,  and  altitude.  There  are 
fourteen  measurements  available  to  the  filter,  including  four  satellite  pscudoranges, 
six  transponder  ranges,  velo  ity  in  3  axes  via  Doppler  aiding  and  altitude  from  the 
barometric  altimeter.  Chapter  III  provides  a  detailed  description  of  these  measure¬ 
ment  sources.  Residuals  are  sent  from  the  Kalman  filter  to  the  FDI  system  and 
corrections  can  be  fed  back  to  the  niter.  The  INS  is  inherently  unstable  in  the  alti¬ 
tude  channel  and  receives  stability  aiding  from  a  barometric  altimeter  (l:p.  83).  A 
sampling  period  of  two  seconds  (sec)  was  chosen  for  all  the  measurements.  Previous 
AFIT  research  (23,  26,  27)  used  update  periods  ranging  from  two  to  ten  seconds. 
The  slower  sampling  periods  were  typically  chosen  to  speed  up  the  simulations  which 
took  several  days  to  run.  This  thesis  used  truth  and  filter  models  with  fewer  states 
resulting  in  shorter  simulations  so  a  higher  sampling  period  was  feasible.  However, 
data  files  were  used  to  store  several  values  at  each  sample  time  and  these  files  became 
unmanageable  with  sample  periods  below  two  seconds.  A  few  simulations  were  also 
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conducted  with  a  one  second  sampling  period  to  ensure  that  no  significant  losses 
occurred  in  the  filter  and  FDI  performance. 


Corrected 


Figure  1.3.  NRS  Block  Diagram 

3.  The  NRS  will  be  modeled  using  differential  equations  that  describe  the  physical 
relationships  between  the  real  world  and  the  NRS  electronics.  The  model  will  be 
based  on  error  states  rather  than  actual  states  because  error  states  will  provide  more 
accurate  estimates  of  position  and  velocity  [2).  Research  at  CIGTF  is  based  on  a 
Litton  LN-93  model  for  the  INS  (9),  a  generic  GPS  model  (23),  and  a  simplified 
RRS  model  (23).  The  models  used  for  this  thesis  are  fully  described  in  Chapter  III. 

4.  Failures  will  be  assumed  additive  rather  than  multiplicative,  allowing  them  to  be 
represented  as  an  additional  term  in  the  equations  describing  the  measurements  from 
the  GPS  rather  than  changes  in  the  dynamics  model.  This  assumption  simplifies  the 
problem  significantly  without  a  loss  of  realism.  Finally,  it  will  be  assumed  that 
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only  one  failure  can  occur  at  a  time  and  that  multiple  failures  are  unmodeled.  See 
Sections  1.6.1  and  3.4. 

5.  Computer  simulations  will  be  run  using  a  program  called  Multimode  Simulation  for 
Optimal  Filter  Evaluation  (MSOFE)  (22).  MSOFE  is  specifically  used  to  provide 
analysis  of  designs  involving  Kalman  filters.  A  complete  analysis,  short  of  actual  test 
flights,  is  possible  using  Monte  Carlo  run  simulations  with  MSOFE.  The  flight  profile 
will  be  generated  using  a  computer  program  called  PROFGEN  (21).  PROFGEN 
takes  user  input  commands  that  describe  aircraft  maneuvers  and  produces  computer- 
compatible  data  which  is  fed  to  MSOFE  for  the  simulations. 

6.  The  state  dynamics  matrix  F  is  considered  piecewise  constant  between  sample  storage 
times  of  two  seconds.  The  discretization  process  in  Matrix*  (8)  takes  this  matrix  and 
uses  an  8th  order  Pade  approximation  for  the  matrix  exponential  EXP(FAt).  See 
Sections  2.2.2  and  4.1. 

7.  The  measurements  are  assumed  sufficiently  uncorreiated  so  the  off-diagonal  terms 
of  the  measurement  noise  matrix,  R,  are  negligible  and  can  be  set  to  zero.  The 
residual  covariance  matrix  A  is  assumed  nearly  diagonal  (or  at  least  dominated  by 
its  diagonal  terms)  based  on  the  relationship: 

A(y)  =  H(/,)P(<DHT(t<)  +  R  (11) 

where  R  is  assumed  stationary  because  it.  has  very  smaii  variations  throughout  the 
simulation.  Although  a  diagonal  R  does  not  ensure  a  diagonal  A,  Equation  (1.1) 
indicates  that  observation  of  the  diagonal  terms  of  A  will  provide  insights  into  the 
condition  of  the  measurements. 

8.  The  simulation  and  post-processing  software  (MSOFE  and  Matrixx)  are  coded  to 
run  in  double  precision  in  order  to  handle  large  disparities  in  various  Kalman  filter 
values. 

9.  MSOFE  simulations  for  filter  tuning  are  performed  using  10-run  Monte  Carlo  anal¬ 
yses  with  statistical  values  averaged  over  the  10  runs.  Single-run  Monte  Carlo  sim- 
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illations  arc  used  when  testing  the  failure  detection  algorithm  and  will  be  compared 
with  multi-run  results  to  ensure  the  single  run  is  not  unusual. 

10.  Failure  thresholds  will  be  determined  empirically  based  on  simulation  results  in  var¬ 
ious  situations.  Sec  Section  2.3.3. 

11.  Taylor  series  truncated  to  first  order  will  be  used  for  linearizing  nonlinear  equations. 
Perturbations  about  some  nominal  point  will  be  established  in  each  case.  See  Sections 
2.2.1,  3. 3. 2.2  and  3.3.3.2. 

12.  The  FDI  algorithm  will  only  view  the  data  within  a  set  window  of  time  in  order  to 
avoid  a  growing  set  of  hypotheses,  as  discussed  later  in  Section  1.5.3.  The  size  ol  the 
window  will  be  determined  empirically  based  on  simulation  results  and  the  window 
will  slide  in  time  to  cover  all  the  data  in  the  simulation.  For  example,  a  ten-sample 
window  would  have  new  data  added  to  the  end  of  the  window  and  simultaneously 
have  old  data  deleted  from  the  beginning  of  the  window.  Additionally,  the  failure 
is  assumed  to  occur  at  the  beginning  of  the  window  to  simplify  calculations.  The 
consequence  of  this  simplification  is  a  delay  in  detecting  the  failure  caused  by  waiting 
for  the  failure  to  reach  the  beginning  of  the  window,  as  discussed  in  Section  2.3.1. 

13.  A  Doppler  system  is  available  to  provide  velocity  aiding  to  the  INS.  The  measure¬ 
ments  from  the  Doppler  are  ideal  and  tell  the  filter  the  exact  error  between  the 
filter  state  and  the  truth  state.  This  ideal  situation  was  assumed  to  allow  direct 
comparison  of  results  against  those  obtained  in  previous  AFIT  theses.  Sec  Section 
3. 3. 1.3. 

1.5  Literature  Review 

This  section  contains  a  review  of  literature  pertaining  to  four  FDI  techniques  with 
specific  interest  on  their  application  to  the  integrated  NRS  described  in  Section  1.4.  More 
details  on  this  NRS  can  be  found  in  a  variety  of  sources,  including  Air  Force  Institute 
of  Technology  (AFIT)  theses  (23,  27).  References  made  to  the  Kalman  filter  and  FDI 
system  will  directly  apply  to  the  components  shown  in  Figure  1.3.  Kalman  filter  theory 
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is  presented  by  Maybcck  (13).  The  literature  review  conclusion  will  discuss  the  validity  of 
each  technique  for  use  in  the  thesis  research. 


1.5.1  Simple  Erroi'  Detection.  A  simple  method  for  error  detection  is  to  monitor 
the  errors  in  the  position  and  velocity  calculations.  A  failure  has  occurred  if  these  errors 
exceed  an  established  failure  threshold,  as  shown  in  Figure  1.2.  The  problem  with  this 
method  is  that  the  desired  values  of  true  position  and  velocity  are  not  known,  so  it  is 
not  possible  to  determine  the  inaccuracy  of  the  calculations.  An  alternative  is  to  use  the 
Kalman  fdter  to  provide  a  statistic  known  as  code  loop  pseudorange,  tracking  error  which 
is  related  to  the  desired  position  through  differential  equations  (24:pp.  1460-1461).  The 
tracking  error  is  sent  to  the  FDI  system  for  comparison  against  a  failure  level.  The  failure 
level  will  be  determined  based  on  past  experience,  flight  tests,  and  most  commonly  through 
computer-generated  Monte  Carlo  simulations  or  covariance  analyses  (ll:pp.  102-106).  If 
the  code  loop  error  tor  a  given  satellite  is  larger  than  the  failure  level,  the  GPS  is  considered 
to  have  faaled  or  lost  lock  m  the  channel  associated  with  that  satellite  and  the  G 1  S  will 
no  longer  use  information  from  this  channel  until  the  satellite  has  regained  lock. 


1.5.2  Direct  and  Analytic  Redundancy.  One  of  the  simplest  and  most  reliable  fail¬ 
ure  detection  techniques  is  the  use  of  redundant  elements  for  voting.  Given  a  system  with 
triple  redundancy,  an  algorithm  can  be  easily  written  that  will  compare  the  outputs  of 
each  element,  allowing  them  to  vote  on  the  condition  of  the  other  elements.  Simply  stated. 


totally  different  value,  the  latter  is  considered  inadequate  to  provide  accurate  information 
and  is  removed  from  the  system.  Once  the  third  element  is  removed  from  the  system,  the 
algorithm  is  unable  t,o  isolate  a  failure.  If  the  two  remaining  elements  disagree,  a  failure 
has  been  detected  but  not  isolated  and  both  elements  will  have  to  be  removed  from  the 
system  or  the  system  performance  will  be  degraded.  The  major  disadvantage  of  direct 
redundancy  is  the  need  for  redundant  hardware  (5:pp.  1-2). 

A  more  sophisticated  approach  uses  analytic  redundancy.  By  using  the  physical  and 
dynamic  relationships  between  instruments  on  an  aircraft,  it  is  possible  to  generate  mul¬ 
tiple  sources  of  the  same  information  mathematically.  These  sources  are  used  by  the  FDI 
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structure  through  voting  techniques.  Analytic  redundancy  avoids  the  need  for  excessive 
hardware,  because  each  instrument  is  required  for  specific  mission  requirements.  The  sur¬ 
plus  of  information  is  also  used  to  provide  isolation  of  the  failure  without  the  need  for 
triple  redundancy  (5:pp,  3-7), 


1.5.3  Multiple  Generalized  Likelihood  Ratio  Testing.  In  some  situations  direct  re¬ 
dundancy  is  not  practical.  An  alternative  is  to  use  a  Kalman  filter  as  shown  in  Figure  1.3 
to  compensate  for  failures  in  the  NRS  and  provide  accurate  estimates  of  the  position  and 
velocity.  Figure  1.4  shows  how  the  sensor  (GPS,  INS,  or  RltS),  the  Kalman  filter,  and  the 
FDI  block  of  Figure  1.3  interact.  Three  hypotheses  are  considered  with  IIo,  Hi,  and  Il2, 
representing  no  failure,  jump  failure  and  ramp  failure  respectively.  The  Kalman  filter  is 
designed  based  on  Hu,  and  two  matching  filters  are  designed  based  on  Hi  and  i/2  (29). 
When  considering  H i  and  //2,  design  parameters  within  the  matching  filters  determine 
what  type  of  failure  is  being  matched,  but  the  magnitude  of  the  jump  and  rate  of  the 
slope  do  not  have  to  be  predetermined  and  arc  estimated  by  the  GLR  algorithm.  1  he 
matching  filters  monitor  the  residuals  provided  by  the  Kalman  filter,  and  each  matching 
filter  computes  a  generalized  likelihood  ratio  (GLR).  The  GLll’s  are  indications  of  which 
hypothesis  is  most  correct.  The  GLR’s  use  maximum  likelihood  estimates  (MLE)  and  are 
compared  through  test  logic  to  detect  and  isolate  failures.  A  corrective  signal  can  be  fed 
back  to  the  Kalman  filter  for  adaptation  to  the  failure.  In  many  cases,  simply  adding  or 
subtracting  a  bias  to  the  sensor  allows  the  system  to  continue  operation  without  losing  the 
sensor  in  question.  One  advantage  of  the  GLR.  test  over  other  FDI  algorithms  is  that  prior 
knowledge  of  the  magnitude  of  the  failure  is  not  necessary.  Further  detail  on  GLR’s  and 
MLE’s  are  presented  by  Willsky  and  Jones  (30,  31). 


The  simple  hypotheses  for  failure  and  no-fail  conditions  may  not  provide  the  ro¬ 
bustness  needed  to  detect  and  adapt  to  certain  failures.  Detection  of  ramp  failures  is 
particularly  difficult  w'heu  the  ramp  rate  varies  significantly.  Figure  1.2  shows  the  time 
delay  involved  in  detecting  ramp  failures.  This  delay'  occurs  because  the  deceiving  signal 
is  slowly  moving  away  from  the  desired  signal  and  takes  more  time  to  cross  the  failure 
threshold.  By  adding  filters  designed  to  match  the  ramp  failure,  these  longer  delays  can 
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Figure  1.4.  Multiple  GLR  Testing 

be  avoided.  It  is  also  desirable  to  design  the  filter  based  on  the  time  of  the  failure.  If  the 
filter  is  customized  to  look  for  a  failure  at  a  specific  instant  in  time,  then  it  has  a  better 
chance  of  good  detection.  This  idea  results  in  a  very  large  bank  of  filters  for  long  periods 
of  time.  The  disadvantage  of  adding  filters  is  the  increase  in  computations  required  for 
multiple  filters.  A  solution  to  this  problem  uses  a  set  number  of  filters  over  a  window  of 
time.  The  number  of  filters  remains  constant,  and  the  performance  of  the  FDI  system  is 
often  maintained  (31:pp.  6U6-6U8).  The  justification  for  using  sliding  windows  is  discussed 
by  Willsky  (31:pp.  G04-605). 

1.5.4  Chi-Square  Testing.  Another  FDI  method  based  on  residual  monitoring  is 
a  chi-square  test  which  is  similar  to  GLR  testing  in  that  it  calculates  a  random  variable 
x(k)  based  on  the  filter  outputs.  One  primary  difference  between  these  two  algorithms  is 
that  GLR  tests  are  functions  of  the  dynamics  model  and  chi-square  tests  are  not,  as  shown 
mathematically  in  Section  2.3.  A  second  major  difference  is  that  clii-square  tests  do  not 
try  to  match  the  failure  and  only  have  one  hypothesis,  making  it  a  binary  test  for  fail/no¬ 
fail.  Without  the  ability  to  distinguish  between  types  of  failures,  the  chi  -square  test  is  not 


MO 


good  for  isolation.  However,  chi-square  testing  is  easy  to  implement,  runs  quickly,  and  can 
provide  the  first  level  of  detection  in  a  multi-level  FDI  scheme. 


1.5.5  Multiple  Model  Adaptive  Estimation.  A  final  FDI  technique  for  discussion  is 
the  use  of  multiple  models  that  represent  the  dynamics  of  the  aircraft  under  a  variety  of 
conditions.  Although  this  technique  is  analogous  to  multiple  GLR  testing  in  many  ways, 
it  differs  in  its  structure  and  decision  making  process.  Unlike  the  multiple  GLR  testing 
which  modeled  a  variety  of  failures  using  matching  filters,  the  multiple  model  adaptive 
estimation  (MMAE)  technique  models  the  dynamic  nature  of  the  aircraft  and  its  sensors 
to  represent  their  behavior  in  the  presence  of  a  failure.  Figure  1.5  shows  that  a  separate 
Kalman  filter  is  designed  for  each  failure  condition,  and  the  residuals  are  used  to  determine 
which  filter  best  models  the  aircraft  and  its  sensors  at  the  current  time.  The  Kalman 
filter  for  the  jump  or  ramp  might  actually  consist  of  many  filters  designed  for  various 
magnitudes  and  slopes.  A  probability  of  accuracy  ranging  from  zero  to  one  is  computed 
lor  eacn  id  ter  anu  multiplied  by  the  filter  estimates  of  position  and  velocity1  to  weight  them 
appropriately.  The  probability-weighted  estimates  are  added  together  1  ■  form  blended 
estimates.  This  blending  allows  for  partial  failures  in  a  sensor  or  combinations  of  failure 
types.  A  probability  of  one  indicates  that  a  filter  is  100%  accurate  in  its  modeling  and  will 
completely  determine  the  final  blended  estimates.  A  probability  of  zero  indicates  that  a 
filter  is  completely  inaccurate  in  its  modeling  and  will  not  affect  the  blended  estimates.  The 
probability  weighting  computation  block  in  Figure  1.5  represents  the  FDI  block  in  Figure 
1.3,  and  the  multiple  Kalman  filters  are  in  place  of  the  single  filter.  MMAE  is  described  in 
detail  by  Brown  and  Hwang  (3)  and  is  extensively  used  in  several  AFIT  theses  concerning 
stochastic  estimation  and  control  (7,  16,  17,  20,  28).  An  important  aspect  presented  in 
these  theses  is  the  use  of  multiple  model  adaptive  control  (MMAC)  for  system  stability 
and  failure  correction. 


1.5.6  Literature  Review  Conclusion,  The  simple  error  detection  presented  in  Sec¬ 
tion  1.5.1  would  normally  be  effective  for  detecting  failures,  but  the  tracking  errors  in  the 
code  loop  pseudorange  caused  by  jamming  or  spooling  are  small  compared  to  other  intrin- 
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Figure  1.5.  Multiple  Model  Adaptive  Estimation 


sic  errors  in  the  GPS.  The  performance  of  the  GPS  would  be  degraded,  but  the  failures 
would  not  be  detectable  with  this  FDI  technique  (4). 

Direct  redundancy  techniques  are  not  feasible  because  multiple  GPS  receivers  on  a 
single  aircraft  are  not  practical.  It  would  also  be  very  expensive  to  support  extra  satellites 
in  space  designed  to  provide  redundancy.  The  concept  of  analytic  redundancy  is  inherent  to 
the  NR.S  design.  The  GPS  receiver,  INS  measurement  unit,  transponders,  and  barometric 
altimeter  are  sensors  which  are  mathematically  related  to  generate  multiple  sources  of  the 
same  information. 

Multiple  GLR  testing  and  MMAE  are  prime  candidates  for  the  FDI  system.  Both 
of  these  techniques  possess  the  versatility  needed  for  the  complex  NRS,  but  the  major 
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difference  is  the  type  of  multiple  filters  required  by  each  technique.  Multiple  GLR  testing 
uses  a  bank  of  matching  filters  that  are  somewhat  less  complicated  than  the  bank  of  Kalman 
filters  used  in  MMAE.  Both  techniques  were  considered  as  final  choices  for  the  FDI  system 
design  and  the  GLR  test  was  chosen.  Additionally,  a  chi-square  test  is  investigated  to  see  if 
it  can  provide  information  in  conjunction  with  the  GLR  algorithm.  This  leads  to  a  multi¬ 
level  FDI  scheme  in  which  the  chi-square  tests  provides  accurate  and  quick  detection,  while 
the  GLR  algorithm  emphasizes  isolation  of  the  failure. 

1.6  Methodology 

The  three  main  steps  of  the  research  approach  are  explained  in  Sections  1.6.1,  1.6.2 
and  1.6.3.  Steps  one  and  two  will  satisfy  the  primary  goal  of  the  thesis,  while  step  three  is 
aimed  at  the  follow  -on  goal  included  as  a  recommendation  for  future  research. 

1.6.1  Preliminary  Studies.  As  mentioned  in  the  literature  review,  multiple  GLR 
testing  has  the  advantage  of  needing  only  one  Kalman  filter,  as  compared  to  M  M  A  E  which 
requires  a  bank  of  such  fdters  in  parallel.  Based  on  assumption  4,  a  simplified  GLR  test 
as  formulated  by  Willsky  (31)  can  be  used  under  the  stipulation  that  the  failures  must  be 
additive  in  nature.  This  assumption  may  not  be  valid  in  computer  simulated  tests,  or  more 
importantly,  in  flight  tests.  Therefore,  it  is  crucial  that  a  preliminary  study  be  performed 
on  the  GLR  technique  to  verify  its  usefulness  in  this  thesis. 

A  basic  satellite  orbit  estimation  problem  can  quickly  test  the  GLR  algorithm.  The 
orbit  problem  is  dearly  defined  by  Maybeck  (14:pp.  46-48)  with  range  measured  in  radius 
units  and  time  in  time  units.  This  study  covers  both,  ramp  and  jump  failures  in  the  range 
measurement  at  different  magnitudes  and  rates,  as  shown  in  Table  1.1  for  this  example.  The 
only  goals  of  the  example  study  were  detection  and  identification  of  a  failure  in  the  range 
estimate  of  the  satellite.  No  effort  was  made  to  adapt  to  any  failures.  The  GLR  algorithm 
performed  favorably,  so  the  multiple  GLR  testing  concept  was  used  rather  than  an  MMAE 
configuration.  No  preliminary  study  was  developed  to  test  the  MMAE  technique. 
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Table  1.1.  Orbit  Runs 


Run 

Failure  Type 

Magnitude  or  Rate 

1 

None 

N/A 

2 

Jump 

0.03  radius  units 

3 

J  ump 

0.5  radius  units 

4 

Ramp 

0.75  radius  units/time  unit 

5 

Ramp 

1.0  radius  units/time  unit 

1.6.2  System  Tests.  With  the  FDI  technique  chosen  in  the  preliminary  study,  ex¬ 
tensive  system  tests  were  performed  on  the  NRS.  The  steps  in  these  tests  include: 

1.  Completely  define  the  models  for  the  GPS,  INS,  and  RRS  based  on  prior  research 
by  Negast,  Snodgrass,  Stacey  (23,  26,  27). 

2.  Refine  the  FORTRAN  code  which  represents  these  models  and  execute  it  with 
MSOFE. 

3.  Write  computer  code  that  will  take  the  outputs  of  the  MSOFE  simulation  and  run 
the  multiple  GLR  test  algorithm. 

4.  Perform  multiple  simulations  accounting  for  various  failures  and  analyze  the  results. 

5.  Modify  the  design  to  get  improved  performance  and  return  to  step  2. 

1.6.3  Adaptive  Techniques.  As  mentioned,  the  follow-on  goal  is  to  determine  cor¬ 
rective  feedback  based  on  the  failure  identification.  This  corrective  feedback  will  improve 
the  ability  of  the  NRS  to  adapt  to  failures  or  uncertainties  in  the  modeling  of  the  NRS.  A 
variety  of  techniques  are  available  and  some  ideas  are  discussed  in  Chapter  V  even  though 
time  constraints  prevented  research  into  this  area. 

1.6.4  Stopping  Criteria.  There  are  three  main  criteriafor  termination  that  coincide 
with  the  three  main  steps  in  the  research  approach. 

1.  In  order  to  prevent  excessive  time  being  spent  on  the  preliminary  study  and  jeop¬ 
ardizing  the  overall  research  effort,  a  suspense  was  placed  on  the  pursuit  of  a  work¬ 
ing  algorithm  for  GLR  testing.  If  an  adequate  failure  identification  algorithm  using 
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GLR  testing  did  not  work,  then  multiple  GLR  tests  would  be  abandoned  and  MMAE 
would  be  pursued  for  the  remainder  of  the  thesis  research.  The  GLR  algorithm  was 
successful  in  the  preliminary  studies  and  was  pursued  throughout  the  research. 

2.  With  one  of  the  FDI  techniques  fully  developed,  several  factors  will  determine  when 
the  system  tests  are  complete. 

(a)  The  design  for  the  basic  integrated  NRS  in  a  benign  environment  must  provide 
accurate  positions  and  velocities  of  the  aircraft.  Specifically,  the  NRS  must 
determine  the  aircraft’s  latitude  and  longitude  with  a  la  accuracy  of  13  feet  (ft) 
horizontal  and  40  ft  vertical.  The  1  a  accuracies  for  velocity  are  0.1  ft/second 
(fps)  in  the  north  and  east  directions,  and  0.4  fps  vertically. 

(b)  The  FDI  algorithm  must  be  able  to  identify  failures  consistently.  This  criterion 
applies  to  all  types  of  failures  and  variations  in  the  magnitudes  and  rates  of 
these  failures.  If  this  goal  were  not  reached  for  certain  failure  types,  then  a, 
determination  would  be  made  either  to  continue  work  on  these  failure  types 
or  to  pursue  adaptive  techniques  for  the  failures  that  have  been  accurately 
identified.  The  decision  was  made  to  focus  on  identification. 

3.  Adaptive  techniques  strive  to  meet  two  goals.  First,  stability  of  the  NRS  must  be 
accomplished.  Stability  implies  that  the  errors  in  the  NRS  estimates  of  position  and 
velocity  are  not  growing  unbounded  as  time  goes  on.  The  second  goal  is  to  achieve 
the  same  accuracies  listed  in  2(a)  above. 

1.7  Overview  of  Thesis 

Chapter  II  presents  the  detailed  theory  used  in  the  research.  Kalman  filter  theory  is 
introduced  with  special  attention  on  discretizing  the  dynamics  of  the  sampled  data  Kalman 
filter.  The  basics  of  the  FDI  algorithms  are  discussed,  including  the  equations  implemented 
for  the  GLR  and  chi-square  tests.  Finally,  the  methods  used  for  selecting  thresholds  are 
presented. 

Chapter  III  fully  describes  the  navigation  system’s  parameters  and  operational  details 
through  an  overall  system  description.  Models  for  ihe  NRS  system  to  include  the  INS,  GPS 
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and  RRS  are  defined  in  detail.  The  failure  models  used  in  the  simulations  and  the  GLR 
algorithm  are  also  introduced. 

Results  of  the  work  done  are  shown  in  Chapter  IV.  The  reduced  order  Kalman  filter 
is  analyzed,  and  a  discussion  of  the  FDI  performance  is  presented.  Chapter  V  summarizes 
the  research  through  conclusions  and  recommendations. 
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11.  Kalman  Filtering  and  Failure  Detection 


2.1  Overview 

This  chapter  presents  the  fundamental  theory  for  application  of  a  Kalman  filter 
to  the  navigation  problem.  Basic  filter  equations  will  be  presented  for  continuous-time 
dynamics  followed  by  a  discretization  process  to  convert  the  filter  for  processing  by  the 
GLR  algorithm.  The  equations  for  the  GLR  test  will  then  show  how  the  failure  detection 
is  accomplished.  Threshold  selection  will  be  discussed  along  with  tradeoffs  in  tuning  the 
filter. 

2.2  The  Extended  Kalman  Filter 

An  extended  Kalman  filter  (EKF)  is  chosen  to  provide  state  estimates  depicting  the 
dynamics  of  the  NRS  components.  The  EKF  allows  for  nonlinear,  time-varying  dynamics 
and  measurement  vectors  as  found  in  this  navigation  problem.  These  vectors  are  linearized 
through  approximation  techniques  about  a  nominal  trajectory  to  form  a  linearized  Kalman 
filter  (LKF).  The  LKF  is  the  basis  for  the  EKF,  which  is  found  by  linearizing  about  the 
updated  state  estimate  rather  than  a  nominal  trajectory. 

2.2.1  The  Sampled  Data  Kalman  Filter.  Let  the  system  model  be  expressed  as  a 
state  equation  of  the  form 

x(t)  =  f[x(t),tj  +  G(t)w(t)  (2.i) 

where  the  state  dynamics  vector  f[x(i),  t ]  is  a  nonlinear  function  of  the  state  vector  x(t) 
aDd  time  t.  Let  the  process  noise  input  matrix  G(r.)  =  I  and  w(t)  be  a  white  Gaussian 
noise  with  mean: 

m  W-E  {w(£)}  =  0  (2.2) 

and  noise  strength  Q(i)  defined  by: 

E  {w(f)wT(t  +  r)}  =  Q(i)«(r)  (2.3) 
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Incorporate  measurements  z (£,)  into  the  filter  at  discrete  times  and  define  them  as  a  non¬ 
linear  function  of  the  state  vector  and  time: 


z (U)  =  K[x(<i),  v(t{) 


(2.4) 


where  v(f,)  is  a  zero-mean  white  Gaussian  noise  of  covariance  R(£,-)  defined  by: 


£{v(  ti)vT(t:)} 


for  ti  =  ij 
for  ^  f  tj 


(2.5) 


and  h[x(f<),  ti]  is  the  nonlinear  observation  vector.  The  LKF  is  based  on  perturbation  states 
about  a  nominal  state  trajectory  xn(Z)  satisfying  xn(t0)  =  xno  and 


x„(t)  =  f[xn(t),f]  (2.6) 

where  f[x„(i),ij  is  shown  in  Equation  (2.1).  The  measurements  are  also  based  on  the 
nominal  states  and  defined  as: 

zn{ti)  =  h[xn{ti),ti]  (2.7) 

The  perturbation  states  are  found  by  subtracting  the  nominal  states  in  Equation  (2.6) 
from  the  original  states  in  Equation  (2.1): 


[x(t)  -  x„(<)j  =  f{x{t),t]  -  f[xn(t),/.]  +  G(/.)w(t)  (2.8) 


Equation  (2.8)  is  approximated  to  first,  order  through  a  truncated  Taylor  series  expansion: 


fx(t)  -  F[t]x(t)]6xn(t)  +  G(«)w(t) 


(2.9) 


where  5x(t)  are  the  perturbation  states.  The  definitions  for  G (t)  and  w(t)  are  unchanged 
and  the  new  linearized  dynamics  matrix  F[(;x„(t)]  is  found  by  taking  partial  derivatives  of 
f[x(t),  t ]  with  respect  to  x(t)  and  evaluated  at  the  nominal  values  for  the  trajectory  x„(f)  : 


F[«;xn(0] 


df)  x(t),t] 


dx 


X=X„(t) 


(2.10) 
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The  discrete-time  measurements  are  similarly  approximated  to  first  order  and  in  the  per¬ 
turbed  form: 


t>z(U)  =  H[/,  ;x(/j)]t'xn(/)  +  v(U) 


(2.11) 


and  the  same  linearization  process  is  used  for  the  measurement  matrix  H[tj;xn (£*)],  result¬ 
ing  in: 


H[i;;x„  (£,■)]  = 


dh[x(U),U] 


dx 


(2.12) 


tx  xn(<o 

The  LKF  in  this  thesis  generates  error  state  estimates  <5x(i)  which  can  he  added  to  the 
nominal  states  to  provide  whole  states  estimates  x(t)  in  the  form: 


x(£)  =  xn(t)  +  6x(t) 


(2.13) 


The  EKF  will  now  be  formed  by  linearizing  about  the  state  estimate  x  rather  than 
the  nominal  trajectory  xn.  The  following  equations  use  the  notation  t / £,■  to  represent 
the  time  history  of  a  given  variable  conditioned  on  the  measurements  taken  through  the 
time  interval  [£i?£j+l).  Also,  represents  the  value  after  propagation  but  prior  to  the 
measurement  update  and  tf  corresponds  to  the  value  after  the  measurement  update.  The 
state  estimates  and  covariance  values  P (£/£,)  are  propagated  from  £<  to  ti+l  by  solving  the 
following  differential  equations: 


x(£/£.)  =  f[x(£/ti),£] 


P  (t/ii)  =  F[£;x(£/£0]P  (t/U)  +  P(f/<,)Fl>;*(*/h)]  +  G(£)Q(£)GT  (t) 


where 


df[x(t),t] 


ox 


*=*(:/!.) 


(2.14) 

(2.15) 

(2.16) 


and  initial  conditions  are  given  by: 


x(ti/li)  -  k{tf)  (2.17) 

P(Wh)  =  P  (<+)  (2.18) 
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The  discrete-time  measurements  are  processed  in  the  EKF  through  update  equations: 


K(U)  =  P(t-)HT[i,;x(<r)]  {H[t.;x(tr)]P(V)HTlt.;x(i-)]  +  R(^)}"1  (2-19) 


x(«t)  =  x(it  )  +  K(«i)  {Zi  -  h[x(«,-  ),<,]} 

p(t+)  =  p(*r)  - 


where 


H(<<)  =  H  fejxftr)]  = 


dh[x(<,-),  ii] 


dx 


x=*(<r> 


(2.20) 

(2.21) 

(2.22) 


and  K(i,-)  is  the  discrete-time  Kalman  filter  gain.  Note  that,  for  the  EKF,  the  measurement 
and  dynamics  vecotrs  are  calculated  about  the  last  state  estimate  x(t~ )  rather  than  the 
nominal  trajectory  used  by  a  simple  linearized  Kalman  filter. 


2.2.2  The  Discrete-Time  Kalman  Filter.  In  order  to  utilize  the  filter  outputs  in 
the  GLR  algorithm,  it  is  necessary  to  discretize  the  state  dynamics  matrix  into  a  state 
transition  matrix  (STM),  «&(^,tj_1).  All  other  quantities  of  interest  such  as  K  and  H  are 
already  in  discrete  form.  The  STM  must  satisfy  the  following  differential  equation  and 
initial  condition  (13): 

<*[*(*,  ti-i)]/dt  =  F(t)*(<,  *<_i)  (2.23) 

*(*<-!,*<- t)  =  I  (2-24) 

Defining  At  =  ti  —  and  solving  with  F  assumed  constant  over  At  (see  assumption  6) 
leads  to: 

i)  =  eFAt  (2.25) 

The  state  equation  can  now  be  written  in  the  discrete  form 


tx(ti)  =  &(ti,ti_1)6x(ti~l)  +  Gd(«,-1)wj(ti_l) 


(2.26) 


where  Gd  and  wd  are  discrete-time  representations  of  G  and  w  defined  earlier. 
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2,3  Failure  Detection 


This  section  will  present  the  theory  behind  GLR  and  chi-square  testing  for  the 
purpose  of  failure  detection.  Given  the  Kalman  filter  developed  in  Sections  ‘2.1  and  2.2,  an 
algorithm  can  be  used  to  observe  changes  in  the  residuals.  If  the  changes  are  significant, 
they  will  represent  failures  in  the  system  by  causing  the  GLR  or  chi  square  value  to  exceed 
a  threshold.  Windowing  will  be  applied  as  discussed  in  Chapter  I. 

2.3.1  GLR  Equations.  The  primary  goal  is  to  define  a  likelihood  function  /(Z,-,0) 
that,  when  compared  to  a  threshold,  will  identify  the  onset  of  a  failure  such  as  jamming 
or  spoofing.  Two  hypotheses  are  established  with  a  Kalman  filter  based  on  H0  (no  failure) 
and  matching  filters  based  on  Hx  (a  failure  has  occurred).  The  Kalman  filter  state  equation 
from  Section  2.2  is  in  the  discrete  form 


<5x°(ti)  =  $(ti,ti„l)8x0(ti-l)  +  Gd(<<_i)wd(ti_1)  (2.27) 

with  discrete  measurements  described  by 

£z°(t.)  =  H(/j)tfx°(tj)  +  v(/.<)  (2.28) 

The  matching  filters  will  not  provide  state  estimation  but  a.re  designed  for  failure  detection 
and  will  have  the  form 


Sx\L)  =  #(/<,  ti..l)6x\ti_l)  +  Gd(/,.1)wd(<i_1)  (2.29) 


and 

tfz1^)  =  H(/.j)^x1(il) -1-  v(t,-)  -f  <!(<*)«((,•,  6>)j/  (2.30) 

where 


d(ti) 
n(ti,0 ) 
v 

e 


failure  vector 
failure  function 
unknown  size  of  the  failure 
unknown  time  of  the  failure 
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Comparison  of  Equations  (2.28)  and  (2.30)  shows  that  a  matching  filter  can  assume  a 
failure  in  the  system  by  modeling  it  as  some  variation  in  the  actual  measurement  beyond 
the  variations  caused  by  dynamics  of  the  system.  Although  the  failure  is  modeled  as  a 
bias  on  the  measurement,  this  model  ca.n  also  represent  changes  in  the  states  caused  by 
real  world  anomalies.  Further  definition  of  the  new  failure  quantities  in  Equation  (2.30) 
will  show  how  the  failures  are  modeled  in  the  matching  filter.  The  failure  vector  d(i,)  is 
r^-by- 1  where  r  is  the  number  of  measurements.  The  l’s  in  the  failure  vector  indicate  which 
measurement  devices  are  assumed  to  be  induced  with  a  failure  and  the  other  elements  of 
d(tj)  are  zero.  The  failure  function  n(t,, 6)  tells  the  matching  filter  where  the  failure  is 
assumed  to  occur  within  the  window  (see  assumption  12)  and  the  form  of  the  failure  such 
as  a  step  function.  This  allows  the  generation  of  a  different  GLR  based  on  different  failure 
times.  For  example,  if  the  failure  is  assumed  to  be  a  unit  step  and  to  occur  3  time  units 
from  the  front  of  the  window,  then  0  =  3  and 


n(ti,9)  -  1  for  6  >  3 
=  0  for  9  <  3 


(2.31) 


Finally,  the  size  v  simply  dictates  the  unknown  magnitude  of  the  failure,  whereas  n(t,,  0) 
and  d(<i)  are  predetermined  design  parameters.  It  is  important  to  note  that,  by  not 
making  v  a  predetermined  constant,  it  will  actually  be  estimated  by  the  GLR.  algorithm  as 
shown  later  and  can  be  used  to  provide  corrective  feedback.  Section  3.4  provides  a  detailed 
discussion  of  how  the  failures  are  modeled  using  these  equations. 

In  general,  the  likelihood  function  or  GLll,  /(£,-,  0 ),  is  based  on  maximum  likelihood 
estimates  of  9  and  //  designated  as  9  and  i >,  When  considering  all  possible  values  of  9  within 
the  window,  the  GLR  with  the  largest  value  indicates  the  presence  and  time  of  the  failure. 
The  derivation  shown  by  Riggins  (25:pp.  112-115)  illustrates  that  in  forming  /(t,-,0),  it 
is  inherently  maximized  over  v  but  obtaining  an  MLE  of  9  is  based  on  the  definition  of 
n(<j,0).  With  7i(tj,0)=one  for  all  t,-,  only  a  single  likelihood  function  will  he  calculated 
and  it  will  not  detect  the  failure  until  it  reaches  the  beginning  of  the  window.  Therefore, 
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the  GLR’s  are  based  on  u  but  not  0.  The  primary  reason  for  this  decision  is  to  reduce 
the  computation  time  that  increases  significantly  when  calculating  several  GLR’s  based  on 
different  values  of  6. 

The  Kalman  filter  residuals  7 (t()  are  shown  in  Equation  (2.20)  and  are  defined  by 

7 (U)  =  Zi~  h[x(*r )>  U]  (2.32) 

and  tl  '1  residuals  can  be  expressed  for  each  hypothesis  as 

:  7 {U)  =  7 °(U)  #i  :  7 (f.)  =  7 °(U)  +  m(<i,6>)//  (2.33) 

For  a  Kalman  filter  successfully  tracking  the  true  states,  7 °(ti)  will  appear  as  zero-mean 
white  Gaussian  noise  of  covariance  A  =  HP"HT  +  R.  With  a  failure  induced  on  the  mea¬ 
surements,  a  signal  of  unknown  magnitude,  m(t,,^)i/,  will  also  be  present  in  the  residuals 
with  v  defined  earlier  and  m(tj,0)  presented  momentarily.  It  is  the  goal  of  the  GLR  algo¬ 
rithm  to  identify  this  signal  by  recognizing  variations  in  the  residuals  from  their  normal 
unfailed  values.  The  GLR  tests  are  particularly  good  at  detecting  jumps  in  the  residuals 
with  the  key  being  how  closely  the  matching  filters  model  the  actual  failures.  Section  1.5.4 
stated  that  the  GLR  algorithm  is  a  function  of  overall  system  behavior  (3>  and  H)  and 
Kalman  filter  gain  K.  This  is  shown  mathematically  in  Equations  (2.34)  and  (2.35)  with 
the  derivation  shown  by  Riggins  (25:pp.  112-115).  The  failure  residual  offset,  m {U,B)  is 
found  through 

m(f;,0)  =  H(t,)y(L-,0)  +  d(U)n(tif0)  (2.34) 

w'here  the  recursive  failure  quantity  y(L+i,$)  is  given  by 

y(L+i,0)  -  $(/.,*<_i)[I  -  K(ti)H(ti)]y(«j,0)  -  (2.35) 

With  the  failure  assumed  to  occur  at  the  beginning  of  the  window  (see  assumption  12), 
Equations  (2.34)  and  (2.35)  can  be  simplified  by  setting  n(t;,0)=one  for  all 

y{U+i,9)  =  $(L,L~i)[I-  K(L)H(L)]y(L,0)  -  ^(L,/.i_1)K(<i)d(<1)  (2.36) 
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m (U)  -•  K(/,)y(/,)  +  d(/.j) 


(2.37) 


The  consequence  of  this  simplification  is  a  delay  in  detecting  a  failure  because  the  failure  is 
not  realized  until  it  reaches  the  beginning  of  the  window.  The  combination  of  the  Kalman 
filter  outputs  and  the  matching  filter  model  will  determine  the  magnitude  of  the  likelihood 
function  defined  as: 

l{u,0)  =  Sr(ti,0)C~1(ti,0)S(ti,0)  (2.38) 


where 


S(U,0)  = 


i- 1 


j= i 


given 


A(<>)  =  H(<;)P(f~  )HT((,)  +  R 

and  the  MLE  of  the  unknown  magnitude  of  the  failure,  v,  is  found  by: 


(2.39) 

(2.40) 

(2.41) 


(2.42) 


The  residual  covariance  A (tj)  and  the  residuals  are  combined  with  Equations  (2.34)  and 
(2.35)  or  Equations  (2.36)  and  (2.37)  to  give  a  linear  combination  of  the  residuals  S(th0) 

ie  irol  f**  ( •}  .  /J)  dofin  m  l?nn  Afinno  ( O  "^Q*\  n  ti  rl  (  O  /I  VltlttlKt  S  /lorieirUI 

<111 VI  Cl*  UVtV.llUUliOiuv  »uuuv-  V^«*|5  1/  J  VIOIUIVM  .11  i_7V|mi,uiv/iio  y  \  ^  j  •  I  v  ll 

rule  based  on  a  threshold,  c,  would  be 


l(U,0)  >  c=>  FAILURE 

1{L,6)  <  e=*  NO  FAILURE 


(2.43) 


2.3.2  Chi-Square  Equations.  A  chi-square  test  will  be  used  and  is  based  on  the 
Kalman  filter  residuals  y(); )  which  are  zero  mean  and  white  with  known  residual  covariance 
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A(tj).  The  chi-square  random  variable  x(4)  is  given  by 

X{h)=  £  7T(^)A“1(<i)7(/j)  (2.44) 

j-k-N+l 

with  N  being  the  size  of  a  sliding  window.  Notice  that  the  system  dynamics  are  not 
included  in  Equation  (2.44)  and  that  only  one  failure  hypothesis  is  available.  This  agrees 
with  the  discussion  in  Section  1.5.4  about  the  simplicity  of  the  chi-square  algorithm.  A 
detection  rule  based  on  an  established  threshold  c  would  be 

X(<*)  >  £=*>  FAILURE 

X{tk)  <  e  =>  NO  FAILURE  (2.45) 

2.3.3  Threshold  Selection  and  Filter  Tuning.  All  thresholds  used  in  the  FI)I  logic 
will  be  determined  empirically.  Three  major  concerns  wiii  be  considered  in  selecting  final 
threshold  values.  First,  the  thresholds  must  be  low  enough  to  prevent  delays  in  detecting 
failures.  Similarly,  they  must  be  low  enough  to  prevent  missed  alarms  caused  by  GLR  or 
X  values  dropping  below  the  threshold  while  an  actual  failure  is  still  present.  Finally,  the 
thresholds  must  be  high  enough  to  prevent  false  alarms  caused  by  variations  in  the  GLR  and 
X  values.  These  variations  may  result  from  aircraft  maneuvers  or  unpredictable  changes 
in  the  random  measurement  noise.  If  the  filter  is  tuned  sufficiently  for  both  good  tracking 
and  enhanced  failure  detection,  then  these  variations  or  noise  floor  will  be  relatively  low. 

When  tuning  the  Kalman  filter,  a  major  tradeoff  must  be  made  to  meet  the  goals 
for  state  estimation  and  failure  detection.  Adjusting  the  process  noise  and  measurement 
ncise  values  to  enhance  FDI  capabilities  may  degrade  state  estimation  and  vice  versa. 
This  can  be  seen  in  Equation  (2.41)  by  realizing  that  reductions  in  R  will  cause  reductions 
in  A,  resulting  in  better  possible  monitoring  of  the  residuals.  The  consequence  is  that 
reductions  in  R  will  cause  an  increase  in  resulting  in  a  conservatively  tuned  filter.  The 
primary  goal  of  the  NRS  filter  is  to  provide  a  navigation  solution  within  the  operational 
specifications  listed  in  Chapter  I.  Once  this  goal  is  met,  the  measurement  noise  for  the  filter 
can  be  reduced  to  allow  better  FDI,  while  the  process  noise  is  increased  to  maintain  good 
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tuning  of  the  filter.  Observation  of  the  state  estimates,  measurement  residuals  and  the 
residual  covariance  A (f*)  will  indicate  a  point  of  diminishing  returns  using  this  technique. 
Once  the  filter  has  been  tuned  via  multi-run  simulations  to  provide  optimum  residual 
monitoring  without  seriously  degrading  state  estimation,  single-run  simulations  will  be 
used  to  evaluate  actual  FDI  performance. 

2.4  Summary 

This  chapter  provided  the  theorectical  basis  for  the  remaining  chapters.  The  actual 
models  for  the  Kalman  filter  and  the  system  are  presented  in  Chapter  III  along  with  the 
details  of  the  matching  filter  design  for  the  GLR  algorithm.  Verification  of  the  discretiza¬ 
tion  process  is  shown  in  Chapter  IV.  Also  included  in  Chapter  IV,  are  the  results  obtained 
by  applying  the  equations  discussed  in  Chapter  II,  including  the  filter  performance,  FDI 
performance  and  tradeoffs  encountered  concerning  threshold  selection  and  filter  tuning. 
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III.  Navigation  and  Failure  Models 


3.1  Overview 

This  chapter  describes  the  models  for  each  of  the  three  navigation  systems  and  the 
failure  models.  An  overall  system  description  is  given,  followed  by  detailed  state  and 
measurement  equations.  Finally,  the  methods  used  to  simulate  the  various  failures  are 
shown  as  they  apply  to  the  theory  in  Chapter  II. 

3.2  Overall  System  Description 

A  brief  reiteration  of  the  basic  elements  in  the  system  is  helpful  for  this  discus¬ 
sion.  The  three  navigation  systems  are  the  GPS,  INS,  and  RRS.  There  are  14  measure¬ 
ments  provided  to  the  Kalman  filter,  including  four  satellite  vehicle  (SV)  pseudoranges,  six 
transponder  ranges,  three-axis  velocity  aiding  from  a  Doppler  system  and  altitude  from 
the  barometric  altimeter.  A  total  of  9'7  error  states  about  nominal  values  make  up  the 
truth  model  to  represent  the  real  world.  A  total  of  15  error  states  are  used  for  the  Kalman 
filter  model. 

A  block  diagram  depicting  the  NR5  truth  and  filter  models  is  shown  in  Figure  3.1. 
The  true  aircraft  position  x  is  generated  by  PROFGEN  and  provided  to  each  navigation 
system.  The  SV  positions  are  generated  by  ORBIT  and  combined  with  the  true  aircraft 
position  to  obtain  pseudoranges  for  use  by  the  GPS.  Each  navigation  system  generates 
perturbations  from  the  true  range  and  the  final  difference  measurements  are  then  formed 
by  subtracting  the  GPS  and  RRS  measured  ranges  from  their  corresponding  INS  calculated 
ranges.  The  EKF  propagates  equations  that  represent  the  NRS  and  uses  the  measurements 
to  update  estimates  of  the  error  states.  Finally,  these  estimates  are  used  to  correct  the 
original  INS  indicated  position. 

3.3  Model  Descriptions 

The  truth  model  consists  of  41  INS  states,  26  RRS  states,  and  30  GPS  states.  The 
filter  model  consists  of  11  INS  states,  two  RRS  states,  and  two  GPS  states.  The  following 
sections  will  provide  details  and  justification  for  these  model  selections. 
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X 


Figure  3.1.  Truth  a,nd  Filter  Model  Block  Diagram 


tj . 0 . 1  uyo  lvioaei.  inis  section  presents  the  truth  and  filter  models  used  for  the 
INS.  The  INS  is  a  strapped-down  wander  aximuth  system  that  senses  aircraft  motion  via 
gyros  arid  accelerometers  and  is  used  as  the  primary  source  for  navigation.  A  93-state 
model  is  presented  with  specific  interest  in  the  41  states  kept  for  the  truth  model.  The 
reduced  order  filter  model  is  then  discussed. 


3.S.1.1  The  93-State  LN-93  Error  Model,  The  41  INS  states  were  extracted 
from  an  original  93-State  INS  model  based  on  the  Litton  LN-93  error  state  model  which 
is  described  using  six  categories  oi  states: 

=  [  bxyT  t>x3r  6x.3t  6x4t  t>x5r  /)x6t  ]t  (3.1) 
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where  foe  is  a  93  X  1  column  vector  and; 


<5x,  represents  the  “general”  error  vector  containing  13  position,  velocity,  attitude, 
and  vertical  channel  errors. 

fix?  consists  of  16  gyro,  accelerometer,  and  baro-altimeter  exponentially  time- 
correlated  errors,  and  “trend”  states.  These  states  are  modeled  as  first  order 
Markov  processes  in  the  truth  (system)  model. 

<5x3  represents  gyro  bias  errors.  These  18  states  are  modeled  as  random  constants 
in  the  truth  model. 

8x4  is  composed  of  the  accelerometer  bias  error  states.  These  22  states  are  modeled 
in  exactly  the  same  manner  as  the  gyro  bias  states. 

8x5  depicts  accelerometer  and  gyro  initial  thermal  transients,  '’'he  6  thermal  tran¬ 
sient  states  are  first  order  Markov  processes  in  the  system  model. 

8xg  models  the  gyro  compliance  errors.  These  18  error  states  are  modeled  as  biases 
in  the  system  model. 


The  truth  model  system  state  space  differential  equation  is  given  by: 
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A  full  desription  of  the  submatrices  for  this  equation  is  given  in  the  Litton  LN -93  manual 
(9). 


S.3.1.2  The  4 1-State  INS  Error  Model.  The  large  number  of  states  in  the 
LN-93  model  results  in  a  truth  model  that  is  cumbersome  to  the  simulation  software. 
Studies  by  Negast  and  others  have  shown  that  41  states  are  sufficient  to  represent  the  INS 


truth  model  accurately  (10,  23).  With  only  41  of  the  original  93  stater,  selected  for  use 
in  this  thesis,  the  submatrices  in  Equation  (3.2)  are  reduced  significant  y.  Appendix  A, 
dfables  A.l  and  A. 2  show  the  41  states  chosen  for  this  thesis  and  Appendix  B  presents  the 
equations  and  state  dynamics  noise  values  used  in  the  sv.bmatrices. 

3.3. 1.3  INS  Measurement  Models  The  two  measurements  associated  with  the 
INS  are  the  barometric  altimeter  aiding  and  the  Doppler  based  velocity  aiding.  As  men¬ 
tioned  previously,  the  altimeter  aiding  is  used  to  overcome  the  instability  inherent  i*  the 
vertical  channel  of  the  INS.  The  altimeter  output  is  modeled  as  the  sum  of  the  true  altitude 
ht,  the  total  error  in  the  barometric  altimeter  6 kg ,  and  a  random  measurement  noise  v. 
Similarly,  the  INS  calculated  altitude  is  the  sum  of  the  true  altitude  and  the  INS  error 
in  vehicle  altitude  above  the  reference  ellipsoid,  6h.  A  difference  measurement  is  used  to 
eliminate  the  unknown  true  altitude,  ht,  resulting  in  Equation  3.3. 

£  _  n  .  n  1  r  r  .  c  f  1  . 

OZ  —  [/if  C/iJ  —  [/if  +  Vlig j  -j-  V 

=  [l]^/i  —  +  V  (3-3) 

A  perfect  Doppler  system  provides  velocity  aiding  to  the  INS  based  on  assumption  13.  The 
Doppler  aiding  could  come  from  the  GPS  or  a  simple  radar  system.  This  measurement 
source  is  not  particularly  significant  for  most  of  the  thesis,  but  it  does  allow  the  filter  to 
generate  better  estimates  of  the  velocity  states.  A  simple  model  is  assumed  for  the  Doppler 
measurement.  All  thiee  channels  (north,  east,  and  up)  are  represented  by  the  difference 
between  tic  truth  state  velocity  error,  hVtn  and  the  filter  state  velocity  error,  SVh  as  shown 
in  Equation  (3.4). 

bzi  —  6Vt.  -  SVi  where  i  =  x,y,z  (3.4) 

Although  this  model  seems  somewhat  unrealistic  in  that  it  provides  the  filter  with  an 
ideal  measurement  for  velocities,  it  does  not  skew  the  performance  of  the  EDI  algorithm 
because  these  measurements  are  not  used  in  the  EDI  calculations.  The  primary  reason  for 
including  the  Doppler  measurements  was  for  comparison  of  the  filter  performance  against 
prior  AEIT  theses  (23,  26,  27). 


3.3.1.J  The  Reduced  11-State  INS  Filter  Model.  The  number  of  states  in  the 
Kalman  filter  is  further  reduced  from  the  41 -state  truth  model.  Consideration  was  given 
to  magnitudes  of  the  states  and  their  estimability  when  deciding  which  states  could  be 
eliminated  for  the  filter.  The  reduced  order  filter  had  to  be  tuned  to  compensate  for  the 
eliminated  states  by  adjusting  dynamics  noise  and  measurement  noise  values.  Table  A. 5 
shows  the  11-states  used  for  the  INS  filter  model  and  Tables  B.10  and  B.ll  show  the  final 
tuning  parameters  used  in  the  filter. 

S.S.2  The  26-State  RRS  Error  Model.  The  Range/ Range-Rate  System  (RRS)  nav¬ 
igation  aiding  system  is  modeled  using  26  states  for  six  ground-based  transponders  (23, 27). 
The  RRS  interrogates  the  transponders  collecting  the  electromagnetic  (EM)  signals  they 
emit.  Ihese  signals  represent  position  information  used  to  calculate  a  navigation  solution 
which  supports  the  INS.  Table  A. 3  shows  the  26  states  and  their  corresponding  state  num¬ 
ber  in  the  overall  97-state  truth  model.  There  are  two  states  common  to  each  transponder 
and  four  which  are  modeled  separately  for  each  transponder. 

3.3.2. 1  RRS  Model  Equations.  The  two  common  states  for  the  transponders 
are  a.  result  of  errors  in  user  hardware.  They  appear  as  bias  terms  and  are  modeled  as 
random  constants.  Their  state  equations  are  given  by: 


where 


Xt,r  =  range  equivalent  of  interrogator  bias 

=  velocity  equivalent  of  interrogator  bias 


The  initial  conditions  for  these  states  were  chosen  to  be  consistent  with  previous  AFIT 
research  (23,  26,  27)  and  are: 


(3.6) 


and 


P  6r,<>t/(^o)  — 


1  ft2  0 

0  10~4  ft2 / sec2 


(3.7) 


The  four  states  which  arc  unique  to  each  transponder  represent  the  error  in  the  surveyed 
position  (x,y,  and  zj  of  the  transponders  and  the  atmospheric  propagation  delay  between 
the  transponder  and  the  receiving  aircraft  or  user.  The  position  errors  are  modeled  as 
random  constants  and  the  atmospheric  error  is  represented  by  a  first  order  Markov  process. 
The  state  equations  for  these  error  sources  are  shown  below  with  i  used  to  represent  the 
various  transponders: 
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The  initial  conditions  for  these  states  were  chosen  to  be  consistent  with  previous  AFIT 
research  (23,  26,  27)  and  are: 

"X-x  ,y  ,z  o)  —  6 
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(3.9) 
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(3.11) 
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3. 3. 2.2  RRS  Measurement  Model.  The  system  description  identified  measure¬ 
ment  sources  which  included  RRS  ranges.  The  RRS  measurements  indicate  the  range  from 
the  transponders  to  the  user  and  Figure  3.2  .shows  the  errors  in  the  true  positions.  This 
measurement  is  expressed  in  Equation  (3.13)  as  the  sum  of  the  true  range,  error  sources, 
and  a  random  measurement  noise  v. 


~  Rt  +  atm  +  bRi  4-  v 


(3.13) 


where 


Rt 

&  Rat  m 

SRh 
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RRS  range  measurement,  from  transponder  to  user 

true  range,  from  transponder  to  user 

range  error  due  to  atmospheric  delay 

error  due  to  equipment  bias 

zero-mean  white  Gaussian  measurement  noise 


The  true  range  Rt  is  not  actually  known,  so  a  difference  measurement,  hz  must  be  obtained 
using  t':e  range,  J R/ws,  from  the  transponder  to  the  user.  This  range  is  not  a  measurement, 
but  is  calculated  by  the  INS  using  Equation  (3.14). 
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where  Xa  and  XT  represent  the  user  and  transponder  position  vectors  in  the  earth  centered 
earth  fixed  (ECEF)  frame  respectively.  Another  way  to  write  Equation  (3.14)  is: 

Rins  =  yj(xu  -  xTy  +  (y„  -  yT)2  +  (zu  -  zT)2  (3.15) 


Figure  3.2  shows  that  X„  and  Xr  ace  not  completely  known  and  have  some  error.  Based  on 
assumption  11  with  perturbations  representing  tne  errors  in  X„  and  XT,  Equation  (3.15) 
can  be  approximated  and  written  in  terms  of  the  true  range  and  a  truncated  first-order 
Talyor  series: 
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(3.16) 


The  solution  for  Rins  is  found  by  substituting  Equation  (3.15)  into  Equation  (3.16)  and 
evaluating  the  partial  derivatives  to  get: 
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(3.17) 


Finally,  the  difference  measurement  is  given  as: 
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3. 3. 2. 3  I'he  Reduced  Two-Slate  RRS  Filter  Model ■  The  same  goal  of  reducing 
the  number  of  states  in  the  filter  was  met  with  the  RRS  model.  Research  at  CIGTF  has 
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Figure  3.2.  RR3  Measurements 


shown  that  retaining  only  the  first  two  states,  which  are  common  to  all  the  transponders, 
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equipment  bias  and  are  represented  as  SRt  in  Equation  (3.18).  Previous  AFIT  researchers 
have  kept  all  26  states  in  the  filter  because  state  reduction  was  not  a  major  goal.  The 
scenario  for  this  thesis  requires  a  filter  of  few  states,  so  an  effort  was  made  to  prove  that 
the  two-state  filter  was  adequate  for  navigation  through  comparison  to  the  26-state  model 
results  obtained  by  Negast  (23).  Filter  tuning  included  increases  in  strengths  of  dynamics 
and  measurement  noises  with  final  values  shown  in  Tables  13.10  and  B.ll. 


8.8.3  The  30-Siatc  GPS  Error  Model.  The  third  and  final  navigation  system  is 
based  on  EM  signals  transmitted  from  orbiting  satellites.  Although  similar  in  concept 
to  the  RRS,  the  GPS  is  modeled  somewhat  differently.  This  model  has  been  developed 
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through  research  at  AFIT  and  many  of  its  fundamental  concepts  are  addressed  in  a  variety 
of  sources  (12,  23,  26,  27).  The  dynamics  and  measurement  equations  for  the  full  30-state 
truth  model  are  presented  followed  by  the  reduced  two-state  fdter  model.  A  tabular  listing 
of  the  30-state  model  is  shown  in  Table  A. 4  and  the  two  fdter  states  chosen  are  listed  in 
Table  A. 5. 


3. 3.3.1  GPS  Model  Equations.  There  are  five  types  of  error  sources  that  are 
modeled.  The  first  two  states  represent  the  errors  in  the  user  clock  and  are  modeled  as 
follows: 

xUclkb 

ZUclkir 

where 

xUdkb  —  range  equivalent  of  user  set  clock  bias 

xuciki r  =  velocity  equivalent  of  user  set  clock  drift 
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(3.19) 


The  initial  state  estimates  and  covariances  for  these  states  were  chosen  to  be  consistent 
with  previous  AFIT  research  (23,  26,  27)  and  are: 
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0 

y  xudkM  j 

u 

(3.20) 


Puclkb,Vclkdr(to)  ~ 


9.0  x  10 14 ft2 
0 


0 

9.0  X  1010/i2/sec2 


(3.21) 


Because  these  error  sources  are  a  function  of  the  user  equipment,  they  are  common 
to  all  the  SV’s.  The  remaining  four  types  of  errors  are  unique  to  each  SV,  based  on 
their  individual  equipment  and  their  position  with  respect  to  the  user.  One  SV  specific 
error  source  is  the  code  loop  error.  Although  the  code  loop  is  part  of  the  user  equipment 
shared  by  all  the  SV’s,  its  error  magnitude  is  relative  to  each  SV.  Next  is  the  atmospheric 
interference  with  the  EM  signals  as  related  to  the  ionospheric  and  tropospheric  delay  in 
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the  signals  propagation.  The  code  loop  error,  tropospheric  delay,  and  ionospheric  delay 
are  all  modeled  as  first  order  Markov  processes  with  time  constants  shown  in  Equation 
(3.22).  All  three  are  driven  by  zero-mean  white  Gaussian  noise  with  strengths  shown  in 
Equation  (3.25). 


The  fourth  error  source  is  due  to  inaccuracies  in  the  clocks  on  board  the  SV’s  and 
the  final  error  source  is  based  on  line-of  sight,  errors  between  the  SV’s  and  the  receiver. 
The  models  for  these  states  are  shown  in  Equation  (3.22)  -  (3.25). 
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(3.22) 


POM(«.)  = 


0.25  ft2 
0 
0 
0 
0 
0 
0 


0 

1.0/i2 

0 

0 

0 

0 

0 


and  means  and  noise  strengths: 


0 

0 

1.0/t2 

0 

0 

0 

0 


0 

0 

0 

25  ft2 
0 
0 
0 


0 

0 

0 

0 

25  ft2 
0 
0 


0 

0 

0 

C 

0 

25ft2 

0 


0 

0 

0 

0 

0 

0 

25  ft 2 


(3.23) 


E  {WGPs{t)}  =  0 


(3.24) 
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E  {wars(tW orsit  +  0} 


0.5  0  0  0  0  0  0 

0  0.004  0  0  0  0  0 


0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


0.004  0  0  0  0 
0 


0 

0 

0 


0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


ft2 /see  -Mr)  (3.25) 


3. 3. 3.2  GPS  Measurement  Model.  The  pscudorange  measurements  between 
the  user  and  the  SV’s  is  shown  in  Figure  3.3  as  Rars.  This  measurement  is  the  sum  of  the 
true  range,  several  error  sources  and  a  random  noise: 


E0ps  —  Et  +  3Rci  -|-  6 Rtrop  +  SRion  +  fiRSc n  +  $RUc u  +  n  (3.26) 

where 


llOP$ 

— 

GPS  pseudorange  measurement,  from  SV  to  user 

Rt 

true  range,  from  SV  to  user 

6Ra 

= 

range  error  due  to  code  loop  error 

^  R-trop 

=2 

range  error  due  to  tropospheric  delay 

*>Rion 

= 

range  error  due  to  ionospheric  delay 

fiRsdk 

range  error  due  to  SV  clock  error 

KB..  .. 

v  ■  *'U  CIK 

range  error  due  to  Ugpr  clock  error 

V 

— 

zero-mean  white  Gaussian  measurement  noise 

Because  Rt  is  not  available  to  the  filter,  a  substitution  will  be  made  to  eliminate  this  term 
through  the  same  techniques  used  to  solve  for  the  RRS  measurement  equation.  First,  the 
satellite  position  vector  Xs  and  the  user  position  vector  X„  are  defined  as: 


1 

e 

*  > 

xu 

xs 

x„  =  ■ 

Vu 

to 

II 

Vs 

.  . 

(3.27) 
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then  the  pseudorange  from  the  user  to  the  satellites  is  calculated  by  the  INS  as: 


(3.28) 


An  equivalent  form  for  Equation  (3.28)  is: 


R. INS  \f^V  "b  {j)u  Vs)  "b  (^y 


(3.29) 


Based  on  assumption  11  with  perturbations  representing  the  errors  in  X„  and  Xs,  Equation 
(3.29)  can  be  approximated  and  written  in  terms  of  the  true  range  and  a  truncated  firs t- 
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order  Taylor  series: 


R 


J  NS 


Rt  + 


dR,NS(Xs,Xu) 


dX,_ 


6XC 


(xs,xu)„ 


+ 


dAfWS(Xs,X„) 


dX, 


6X„ 


(Xs  )norr 


(3.30) 


The  solution  for  Rjns  is  found  by  substituting  Equation  (3.29)  into  Equation  (3.30)  and 
evaluating  the  partial  derivatives  to  get: 


-Rfmc  —  Rt  ~~ 


Xs  T'u  1  6x, 


l  I*, 


[ 


ys  -  y,j 


\R, 


hu  - 


X  ~  V 

■  \R,ns\  . 


■  6z,, 


i  Xs-  ...  Xu  .fjX  _l  Vu  c  I  ZS  Zu  ]  c 

\R,Ns\  S  \Rims\  Js  L|Hf„U  * 

Finally,  the  GPS  pseudorange  difference  measurement  is  given  as: 


(3.31) 


bz  —  Rins  Rqps 


xs  ~  xu 

■bxv  - 

'vs  -yu 

•SVu  - 

'zs~Zu 

L  \R:,<S\ 

’**  ~XU 

■6xs  + 

rys  -  Vu ' 

■6ys  + 

* S  Z{J 

.  I^nvsl  . 

!*,„s| 

~  [l]<5/icZ  -  [l]t>Rtrop  ~  [1] 

-  [1]<5 Rsdk  -  [l]^„clt  -  v  (3.32) 


3. 3. 3. 3  The  Reduced  Two-State  GPS  Filter  Model.  Various  research  efforts 
have  shown  that  two  states  provide  a  sufficient  model  for  a  GPS  (10,  23)  The  primary 
argument  is  that  the  errors  modeled  by  the  other  28  states  are  sn  )  when  compared  to 
states  one  and  two  which  are  common  to  all  SV’s.  By  adding  measurement  noise  (increasing 
R.)  and  retuning  the  filter,  the  overall  performance  of  the  NR.S  can  be  maintained  with 
a  significantly  reduced  order  model.  The  final  noise  values  are  shown  in  Tables  13.10  and 
B.ll. 
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3.4  Failure  Models 

This  section  discusses  the  methods  used  to  model  failures  in  the  MSOFE  simulations 
and  the  corresponding  models  used  by  the  FDI  algorithm  for  detection. 

3J,.l  Simulation  Failure  Models.  The  10  different  types  of  failures  that  were  mod¬ 
eled  axe  presented  with  actual  values  in  Table  3.1.  A  short  description  of  each  failure  is  as 
follows: 

1.  Jamming  -  Jamming  is  modeled  as  a  sudden  increase  or  jump  in  the  measurement 
noise  associated  with  all  four  SV’s.  Ihis  failure  is  induced  in  all  SV  signals  because 
the  jamming  is  assumed  to  occur  at  the  receiver,  which  will  affect  all  four  channels 
simultaneously.  The  jamming  noise,  RS,  is  added  to  the  system  model  measurements 
only,  with  values  shown  in  Table  3,1.  These  values  represent  various  levels  of  jamming 
which  result  in  lower  carrier-to -noise  ratios,  C/N0,  of  the  GPS  signal  and  were 
calculated  using: 

RS  _  KxBn  K-Ju Bn 

A2  "  (C/N0)  +  0 C/No)2  {  ’ 

where 

RS  =  jamming  noise  added  to  system  measurement  noise 
A  =  code  modulation  chip  width 
C /N0  =  carrier-to-noise  ratio 
Bn  ~  code  tracking  loop  noise  bandwidth 
F>i  =  one-sided  IF  bandwidth 
Ki  —  code  mechanization  parameter  constant 
K 2  =  code  mechanization  parameter  constant 
The  theory  supporting  Equation  (3.33)  can  be  found  in  an  article  by  Martin  (12). 
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Table  3.1.  Failure  Types  and  Models 


|  Run 

Goal 

. 

Method 

Time  Frame 

Comments 

m 

Baseline 

No  Failure 

N/A 

N/A 

■ 

Heavy 

Jamming 

Jump  in 

Measurement  Noise 

Increase  R 
from  2  to  4000 

2000-4000  sec 

C/N0  =  15  dB-Hz 

2 

Medium 

Jamming 

Jump  in 

Measurement  Noise 

Increase  R 
from  2  to  500 

2000-4000  sec 

C/N0  =  20  dB-Hz 

3 

Light 

Jamming 

Jump  in 

Measurement  Noise 

Increase  R 
from  2  to  30 

2000-  4000  sec 

C/N0  ■■=  25  dB-Hz 

4 

Spoofing 

Bias  on  ONE 
Measurement 

Add  bias=7000 
to  SV1 

2000-4000  sec 

Causes  1  mi 
position  error 

5 

Spoofing 

Bias  on  ONE 
Measurement 

Add  bias=700 
to  SV1 

2000-4000  sec 

Causes  500  ft 
position  error 

6 

Spoofing 

Bias  on  ONE 
Measurement 

Add  bias=700 
to  SV1 

4000-6000  sec 

Causes  500  ft 
position  error 

7 

Spoofing 

Bias  on  ALL 
Measurements 

Add  bias=700 
to  SV1-SV4 

2000-  4000  sec 

Causes  small 
position  error 

8 

Spoofing 

Ramp  on  ONE 
Mcas  rement 

Add  ramp— 2T 
to  SV 1 

1500-6000  sec 

Causes  11/2  miles 
position  error 

9 

Spoofing 

Ramp  on  ONE 
Measurement 

Add  ramp=lT 
to  SV1 

2000  -3000  sec 

Causes  700  ft 
position  error 

10 

GPS  Fail 

Loss  of  ONE 
Measurement 

SV1-0  over 
time  frame 

2000-2200  sec 

Causes  complete 
disruption 
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2.  Spoofing 


(a)  Model  A  -  A  bias  is  added  to  the  measurements  associated  with  one  or  all  four 
SV’s  designated  SV1-SV4.  This  bias  will  be  added  suddenly  in  time  and  the 
values  chosen  represent  a  variety  of  spoofers.  An  intelligent  spoofer  will  have  an 
accurate  estimation  of  the  measurements  being  received  by  the  aircraft  from  the 
satellites  and  can  discretely  add  an  undetectable  bias  to  these  measurements.  A 
less  sophisticated  spoofer  would  have  to  use  a  larger  bias  to  ensure  effectiveness 
in  corrupting  the  measurements  while  running  the  risk  of  being  detected.  A 
simple  method  was  used  to  determine  the  amount  of  bias  added  to  the  pseu¬ 
dorange  and  is  illustrated  through  an  example.  Assume  the  spoofer  wanted 
to  pull  the  aircraft  approximately  one  mile  off  course  in  terms  of  latitude  and 
longitude.  The  net  distance  error  would  be  the  magnitude  of  the  change  in  the 
filter  computed  latitude  and  longitude  error  with  and  without  the  spoofer: 

Distance  error  =  \J(60X  -  f>Ox:,p00jY  +  (S6y  -  S6y:,pooJy  (3.34) 


where 


Mx-.'pooj  =  latitude  error  with  spoofing 
SO,. r.,pooj  =  longitude  error  with  spoofing 

The  altitude  channel  is  not  used  to  compute  the  distance  error  because  altitude 
information  is  readily  available  from  other  instruments  and  the  goal  is  to  draw 
the  aircraft  off  target.  Figure  C.3  shows  the  increase  in  the  filter-computed 
errors  in  all  three  channels  and  the  resulting  distance  error  of  approximately 
one  mile  over  the  time  frame  of  the  failure  for  a  bias  of  7000  ft  on  SVl.  Only 
one  of  the  simulation  runs  will  induce  the  bias  in  more  than  one  SV  signal.  It 
is  considered  unrealistic  that  a  spoofer  would  be  able  to  identify  all  four  SV’s 
selected  by  the  GPS  receiver  consistently  because  the  SV’s  are  chosen  based  on 
their  geographic  relation  to  the  receiver.  The  spoofer  would  not  only  be  more 
complex,  but  would  require  four  receiver/transmitters  to  accomplish  this  task. 
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Finally,  simultaneously  adding  a  bias  to  all  four  SV  measurements  would  appear 
as  an  increase  in  the  user  clock  bias  error  and  would  be  quickly  compensated 
by  the  fiiter.  In  contrast,  a  bias  added  to  only  one  of  the  four  GPS  channels 
should  cause  errors  in  the  navigation  solution  and  degrade  filter  performance 
throughout  the  time  frame  of  the  failure.  A  single  run  will  be  performed  to 
verify  these  predictions  and  show  the  distinction  of  these  failure  types. 

(b)  Model  B  -  An  even  more  subtle  failure  would  be  a  ramped  value  added  to  the 
SV  measurements.  The  ramp  rates  are  somewhat  arbitrary  but  with  the  same 
basic  idea  of  slowly  drawing  the  aircraft  off  course.  Slopes  of  one  ft/sec  and  t  wo 
ft/sec  are  simulated  with  net  position  errors  of  700  ft  after  1000  sec  and  1  1/2 
miies  after  4500  sec,  respectively.  Simulations  were  not  performed  with  a.  ramp 
value  added  to  all  the  measurements  as  explained  above. 

3.  GPS  failure  -  Represented  as  a  loss  of  one  or  all  SV  measurements  and  modeled  by 
setting  the  measurement  z=0.  The  time  frame  illustrates  that  the  SV  will  be  lost 
for  a  period  of  time  but  may  be  reacquired.  Simulations  with  all  four  SV  signals 
removed  will  not  be  done,  as  failing  one  channel  has  the  same  effect  of  completely 
disrupting  the  filter. 

3.4-2  FDI  Failure  Models.  The  multiple  GLR  algorithm  allows  for  various  hypothe¬ 
ses  based  on  the  types  of  failures  being  detected.  A  different  failure  model  can  be  formed 
for  each  matching  filter  to  optimize  detection  and  isolation  of  failures  through  observa¬ 
tion  of  the  MLE’s  and  GLR’s  associated  with  each  matching  filter.  In.  simpler  terms,  the 
matching  filter  with  the  largest  GLR  is  deemed  the  most  correct  in  finding  the  failure. 
Given  a  group  of  different  failures,  a  matching  filter  could  be  generated  based  on  each 
failure  type  or  a  less  complex  FDI  scheme  would  be  to  model  only  the  most  likely  fail¬ 
ures.  By  selectively  choosing  a  bank  of  matching  filters  to  cover  most  of  the  failures,  the 
computational  burden  on  the  computer  running  the  FDI  algorithm  can  be  reduced  from 
one  running  an  excessive  number  of  filters  in  parallel.  This  concept  is  used  in  conjunction 
with  a  simple  chi-square  test  to  form  a  two  level  FDI  scheme.  The  first  level  consists  of 
both  the  chi  -square  test  and  five  matching  filters  to  provide  initial  detection  of  a  failure. 
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The  second  level  is  used  to  isolate  the  failure  based  on  the  five  GLR  outputs  alone.  This 
scheme  benefits  from  the  simplicity  and  reliability  of  the  chi-square  test  and  the  flexibility 
of  the  multiple  GLR.  tests. 

The  failure  quantities  in  Equation  (2.30)  are  specifically  defined  to  represent  the 
matching  filter  models.  The  function  n(i,-,0)  is  equal  to  one  for  all  U,  which  implies  that 
the  failure  occurs  at  the  beginning  of  the  window  (see  assumption  12).  The  failure  vector 
d(L)  is  defined  in  general  as 

SKI 
SV  2 
SK3 
SV4 
TRl 
TR2 
TR‘3 
TRA 
TR5 
TRG  J 

where  SV1  -  SV4  and  TRl  -  TR6  represent  the  satellite  vehicles  and  transponders  re¬ 
spectively  and  will  have  a  value  of  zero  or  one  for  this  simplified  model.  Note  that  the 
Doppler  and  altimeter  measurements  are  not  included  in  these  models  because  their  mea¬ 
surement  residuals  are  virtually  unaffected  by  the  failures  in  question  and  they  slow  down 
the  algorithm  when  included  in  these  equations.  The  plots  of  the  GLR  test  using  all 
14  measurements  (Doppler  and  altimeter  included)  are  not  presented  because  they  look 
essentially  the  same  as  the  GLR  test  using  only  10  measurements. 

The  first  four  matching  filters  assume  a  failure  in  only  one  of  the  four  SV’s  and  tire 
last  filter  models  a  failure  in  all  SV’s  simultaneously,  resulting  in  failure  vectors  of  the 
form: 
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SV1  SV2  SV3  SV4  SVl-4 


d(U)  = 


1 

0 

0 

'o' 

1 

0 

1 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

0 

1 

1 

0 

0 

0 

0 

or 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(3.36) 


These  models  should  perform  well  on  bias  type  failures  and  encounter  some  delay  in  de¬ 
tecting  ramp  type  failures.  This  is  expected  since  the  matching  fdters  are  designed  for 
biases  and  not  ramps,  but  the  ramp  results  will  show  how  well  the  filters  perform  on  other 
types  of  failures.  Further  studies  could  include  more  complex  models  specifically  designed 
to  detect  ramp  failures. 


3.5  Summary 

This  chapter  presented  the  details  for  both  the  navigation  filter  and  failure  models. 
The  basis  for  the  measurement  models  was  discussed  to  help  describe  the  intricacies  of 
the  NRS  design.  The  state  and  dynamics  model  descriptions  illustrate  the  high  degree 
of  nonlinearity  and  time-variance  of  the  system.  The  reduced  order  filter  models  were 
presented,  including  consideration  of  the  critical  job  of  tuning  the  NRS  filter.  The  methods 
used  to  induce  the  failures  in  the  simulations  were  shown  along  with  the  models  for  the 
matching  filters  designed  to  detect  and  isolate  these  failures.  Results  and  analysis  of  these 
simulations  are  presented  in  the  next  chapter. 
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IV.  Result:-:  and  Analysis 


4-1  Assumption  Verification 

One  of  the  primary  concerns  in  this  research  has  been  to  verify  assumption  6  pre¬ 
sented  in  Chapter  I.  The  GLR  algorithm  required  a  discretized  version  of  the  state  dynamics 
matrix  F  which  was  assumed  piecewise  constant  over  the  two  -second  sample  period.  All 
43  time  varying  elements  of  F  were  carefully  scrutinized  to  ensure  that  they  did  not  change 
significantly  between  samples.  Figure  C.2(a)  shows  the  worst  case  element  of  F  in  terms 
of  time  variance  over  the  entire  simulation.  From  this  plot  it  is  hard  to  see  if  the  high 
dynamics  are  being  preserved  over  the  sample  period.  Figure  C.2(b)  provides  a  closer  look 
at  the  changes  occurring  in  this  plot  during  the  most  dynamic  time  period.  The  time  scale 
in  figure  C.2(b)  is  reduced  so  that  values  at  every  two-second  sample  are  clearly  seen. 
From  this  plot  it  is  obvious  that  the  quickly  changing  nature  of  the  F  matrix  is  preserved 
with  a  sample  time  of  two  seconds,  so  the  discretized  STM,  4>,  based  on  assumption  6  will 
be  valid. 

4-2  NRS  Filter  Performance 

A  variety  of  simulations  were  run  on  MSOFE  to  determine  the  final  tuning  values 
chosen  for  the  filter.  Observation  of  state  statistics  and  measurement  residuals  gave  insights 
into  adjustments  of  the  tuning  values  necessary  to  provide  good  state  estimation  while 
enhancing  FDI  potential  as  discussed  in  Section  2.3.3.  The  state  plots  in  Appendix  D 
represent  a  well  tuned  filter  with  the  primary  concern  being  to  reduce  the  mean  error  for 
each  state.  A  quick  method  for  evaluating  how  weil  the  filter  is  tuned  is  to  ensure  that  the 
meanisigma  for  the  state  error  is  bounded  by  the  fdter-computed  zeroisigma.  The  FDI 
potential  is  evaluated  using  the  residual  plots  and  the  goal  is  to  bound  the  residualisigma 
by  the  filter-computed  zeroisigma  derived  from  the  residual  covariance.  The  baseline 
plots  shown  in  Figures  D.16-D.20  indicate  that  further  reduction  of  the  measurement 
noise  covariance,  R,  would  improve  the  FDI  performance  based  on  A  =  HP~HT  -|-  R. 
However,  lower  values  of  R  caused  serious  degradation  of  filter  tuning  for  state  estimation 
performance.  It  is  also  unrealistic  for  R  to  have  values  below  the  actual  measurement 
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device  noise  variance  values.  Therefore,  a  tradeoff  was  made  resulting  in  final  values  for  Q 
and  R  shown  in  Tables  B.10  and  B.ll. 

The  15-state  filter  performance  is  compared  against  the  results  of  the  thesis  by  Negast 
(23)  and  the  operational  specifications  listed  in  Section  1.6.4.  All  results  are  based  on  10- 
run  Monte  Carlo  simulations  with  ensemble  averaging  performed  over  the  10  runs,  Negast 
conducted  studies  on  a  variety  of  reduced  order  filters  with  some  of  his  results  shown  in 
Table  4.1.  The  values  presented  are  the  temporal  averages  of  the  ensemble  averages  of  true 
filter  estimation  errors  (la)  for  the  position,  velocity,  and  attitude  errors  over  the  two-hour 
flight  profile.  Both  the  97-  and  69-state  filters  were  evaluated  against  a  128-state  truth 
model  defined  by  Negast.  The  15-state  filter  used  in  this  thesis  was  evaluated  against 
the  97-state  truth  model  defined  earlier.  Although  this  disparity  could  slightly  skew  the 
results,  a  comparison  of  the  filter  performance  is  still  worthwhile  in  demonstrating  the 
ability  of  the  15-state  filter  to  provide  a  good  navigation  solution. 


Table  4.1.  Temporal  Averages  of  True  Filter  Errors  (la) 
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All  the  operational  specifications  were  easily  met  and  exceeded  as  shown  in  Table 
4.1.  It  is  clear  from  the  results  that  the  overall  performance  of  the  15-state  NRS  filter  is 
degraded  in  comparison  to  the  higher  order  filters.  Some  of  the  errors  are  actually  better 
for  the  15-state  filter  than  for  the  69-state  filter  as  a  result  of  either  the  difference  in 
truth  models  or  good  filter  tuning.  Investigation  into  the  individual  state  plots  shown  in 
Appendix  D  reveals  other  problems  encountered  with  the  15-state  filter.  Large  variations 
in  the  true  system  behavior  such  as  high  dynamic  maneuvers  will  test  the  robustness  of 
the  filter.  With  a  reduced  order  filter  like  the  15-state  NRS,  a  single  sharp  turn  or  altitude 
change  is  easily  handled,  but  several  of  these  maneuvers  caused  the  filter  estimation  error 
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to  increase  significantly.  The  flight  profile  shown  in  Figure  C.l  indicates  that  problems 
may  occur  between  4800  and  6600  seconds.  Analysis  of  the  filter  performance  verified 
that  two  or  three  high-g  turns  in  rapid  succession  provided  justification  for  increasing  the 
tuning  values  of  certain  states  in  order  to  enhance  tracking.  Adaptive  tuning  techniques 
were  employed  on  states  one,  two,  12  and  14  with  upper  and  lower  process  noise  values 
shown  in  Table  B.10.  The  upper  noise  values  would  be  employed  following  a  series  of  quick 
turns  and  the  lower  noise  values  would  be  used  again  a  few  hundred  seconds  after  the  last 
major  turn.  A  study  would  have  to  be  done  with  a  variety  of  flight  profiles  to  determine 
the  optimal  conditions  for  selecting  the  upper  or  lower  tuning  values.  The  key  issue  is 
that  adaptive  tuning  is  required  and  was  implemented  for  this  simplified  filter  model  as 
shown  in  Figures  D.l,  D.2,  D.12  and  D.14.  Statistics  for  the  other  11  states  are  plotted  in 
Appendix  D  and  indicate  that  all  the  states  were  estimable  at  various  degrees  of  tuning. 
State  3  represents  the  azimuth  error  and,  as  anticipated,  was  difficult  to  estimate  in  this 
wander-azimuth  system.  Aiso,  its  magnitude  is  very  small  resulting  in  little  or  no  impact 
on  the  NRS  filter  performance  so  fine  tuning  was  not  attempted. 

Another  goal  was  to  verify  the  validity  of  the  two-state  RRS  filter  model.  Although 
temporal  averages  were  not  presented  by  Negast  for  his  26-state  RRS  idter,  visual  exami¬ 
nation  of  his  plots  reveals  RRS  range  bias  errors  of  0.7  to  1.0  ft  in  comparison  to  8.27  ft  of 
error  with  the  two-state  model.  This  seemingly  large  increase  had  a  small  effect,  on  overall 
filter  performance  and  the  two-state  model  is  deemed  to  be  sufficient  for  most  applications. 
The  major  impact  is  the  reduction  in  computational  loading  realized  by  removing  24  states 
from  the  filter. 

4-3  FDI  Performance 

Both  the  GLR  and  chi-square  FDI  algorithms  are  discussed  with  run  numbers  corre¬ 
lating  to  Table  3.1  Specific  concern  is  given  for  missed  and  false  alarms  and  methods  for 
enhancing  failure  detection  are  discussed.  Three  types  of  plots  are  discussed  to  illustrate 
the  performance  of  the  FDI  algorithms.  First,  likelihood  function  plots  show  the  actual 
results  of  the  GLR  algorithm.  Second,  a  fail  flag  routine  was  written  to  smooth  the  GLR 
plots  based  on  a  detection  threshold  and  a  parameter  called  number Jow.  This  routine 
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assumes  there  is  a  failure  ( fail  flag— one)  until  the  GLR  value  falls  below  the  threshold  for 
a  number  of  successive  samples  >  number -low,  resulting  in  fail  jlag=Q.  The  consequence 
of  a  large  number Jow  is  that  fail  flag— one  for  several  samples  even  after  the  failure  has 
been  removed.  Chi-square  plots  are  also  generated  for  three  different  window  sizes.  Tech¬ 
niques  employing  two  fail  flag  routines  with  a  large  and  small  number  Jow  could  help  make 
quicker  decisions  in  dropping  the  fail  flag.  Other  decision  logic  might  be  used  to  overcome 
problems  in  detecting  certain  types  of  failures,  but  this  thesis  focused  on  a  simple  scheme 
using  a  single  fail  flag  routine  and  a  single  chi-square  test  for  decisions. 

4-3.1  Jamming  Detection.  Three  levels  of  jamming  were  induced  as  indicated  by 
runs  one,  two  and  three  in  Table  3.1.  The  degradation  in  the  filter  performance  is  di¬ 
rectly  related  to  the  amount  of  jamming  noise  induced.  Selected  state  plots  are  shown  in 
Appendix  E  to  illustrate  the  effects  of  heavy  jamming  on  the  filter.  The  filter  is  able  to 
reacquire  a  good  estimation  of  the  states  quickly  following  removal  of  the  jamming  noise, 
demonstrating  the  ability  of  the  filter  to  survive  a  hostile  environment.  The  residual  plots 
shown  in  Figures  E.4  and  E.5  make  it  clear  that  simple  residual  monitoring  would  be  suf¬ 
ficient  to  detect  the  jamming  when  compared  with  the  baseline  residuals  in  Figures  D.16 
-  D.17.  Figure  E.6  seems  to  indicate  that  the  transponders  are  not  significantly  affected 
by  the  jamming,  but  even  these  minor  variations  from  the  baseline  in  Figure  I). 18  will  be 
detectable  by  both  the  GLR  and  chi-square  algorithms. 

The  results  of  the  GLR  test  for  heavy  jamming  are  shown  in  Figures  E.7  and  E.8 
and  should  be  compared  to  the  baseline  plot  in  Figure  I). 21.  The  five  plots  in  each  figure 
represent  the  GLR’s  based  on  the  five  different  matching  filters  discussed  in  Section  3.4.2. 
As  expected,  the  matching  filter  assuming  a  failure  in  all  four  satellites  provides  the  best 
detection.  Figure  E.8  shows  that  with  a  threshold=120  and  a  number  Jow=23,  the  fail  flag 
remains  up  (equal  to  one)  throughout  the  failure  time  and  missed  alarms  can  be  avoided. 
The  relatively  large  value  for  number  Jow  is  attributed  to  the  continuous  variations  in  the 
residuals  resulting  from  the  addition  of  random  noise.  These  variations  prevent  the  GLR 
algorithm  from  reaching  a  constant,  level  above  the  noise  floor  and  the  large  number  Jow 
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causes  an  excessive  46-second  delay,  which  is  excessive,  in  determining  that  the  jamming 
is  removed. 

Comparison  of  Figures  D.21  and  E.9  show  that  the  chi-square  test  was  very  successful 
in  detecting  the  jamming.  The  best  choice  for  the  window  size  is  15  samples  in  order  to 
reduce  delays  in  detection  and  false  alarms  which  are  both  apparent  at  the  end  of  the 
failure  time.  A  larger  window  size  would  only  be  necessary  if  the  threshold  had  to  be  set 
so  low  (for  good  detection)  that  it  caused  false  alarms  due  to  the  noise  floor. 

Similar  results  were  obtained  with  the  medium  jamming  level,  as  shown  in  Figure 

E. 10.  Both  the  GLR  and  x  algorithms  are  capable  of  detecting  this  level  of  jamming. 
Neither  algorithm  was  successful  at  consistently  detecting  the  light  jamming,  as  shown  in 
Figure  E.ll.  The  main  reason  the  GLR  test  had  difficulties  detecting  jamming  failures  is 
that  the  matching  filters  are  designed  to  isolate  a  bias  failure  and  the  jamming  noise  is 
quite  different  from  a  bias.  In  contrast,  the  chi-square  test  is  not  dependent  on  a  failure 
model  and  is  better  suited  to  detect  the  variations  in  the  residuals  caused  by  jamming. 

4-3.2  Bias  Failures.  Runs  four  -  six  induced  bias  failures  on  SV1  only.  Figures 

F. l  and  F.'2  show  the  degradation  in  the  filter  throughout  the  time  frame  of  the  failure 
resulting  from  a  bias  of  7000  ft  being  added  to  the  pscu.lorange  of  SV 1.  The  GPS  residuals 
shown  in  Figure  F.3  clearly  indicates  that  a  failure  begins  at  2000  sec  and  falls  ofF  at  4000 
sec.  Although  somewhat  masked  by  the  scale  of  the  plot,  Figure  F.3(b)  shows  that  the 
residual  values  are  also  large  throughout  the  entire  failure  time.  This  is  easier  to  see  on  the 
transponder  residuals  in  Figure  F.4.  An  interesting  observation  is  the  significant  variation 
in  the  transponder  residuals  even  though  no  failure  was  induced  on  their  measurements. 
Similarly,  all  the  SV  residuals  are  affected  as  a  result  of  a  failure  on  SVl  only.  This 
apparent  coupling  of  the  measurements  is  not  due  to  correlation  in  the  real  world,  but 
rather  a  result  caused  by  the  design  of  the  NRS  filter.  The  measurement  models  are 
functions  of  the  error  states  within  the  filter  and  with  large  variations  in  these  estimated 
states,  all  the  measurement  residuals  are  affected.  This  coupling  of  the  measurements 
will  have  a  negative  effect  on  failure  isolation.  It  is  also  important  to  recognize  the  large 
spike  occurring  at  approximately  2300  sec,  making  the  presence  of  the  failure  seem  obvious. 
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Thfs  was  determined  to  result  from  a  sharp  turn  by  the  aircraft  which  dithered  the  system, 
indicating  that  FDI  decisions  should  coincide  in  time  with  dithering.  The  significance  of 
this  dither  signal  will  be  discussed  later. 

A  typical  GLR  plot  is  shown  in  Figure  F.5.  The  lower  plot  has  been  capped  off  at  200 
to  allow  viewing  of  the  data  close  to  the  noise  floor.  This  noise  floor  is  the  result  of  initial 
transients,  variations  in  random  noise  and  changes  in  the  aircraft  dynamics.  A  threshold 
could  be  set  for  quick  detection  and  to  prevent  any  false  alarms  (recognized  as  the  GLR 
value  crossing  the  threshold  when  no  failure  is  present),  but  a  missed  alarm  would  result. 
The  approach  would  be  to  run  the  fail  flag  routine  with  number Jow— two,  resulting  in  a 
two  -second  delay  in  realizing  that  the  failure  is  gone,  no  missed  alarms  and  a  detection 
time  of  two  seconds.  This  concept  is  applied  repeatedly  throughout  the  rest  of  the.  analysis. 

The  GLR  values  for  all  five  matching  fdters  are  shown  in  Figure  F.6.  Analysis  of 
this  data  generates  fail  flag  results  shown  in  Figure  F.7.  It  is  reassuring  that  the  best 
results  are  obtained  from  the  matching  filter  SVl  which  assumes  a  bias  failure  in  SVl. 
There  is  concern  for  the  false  alarms  in  SV2  and  SV4.  One  solution  would  be  to  raise  the 
threshold  above  the  noise  floor  of  both  matching  filters.  Although  no  increase  in  detection 
delay  will  occur  for  this  failure,  later  simulations  with  rarnp  failures  will  suffer  fiom  raising 
the  threshold  too  high  and  causing  delays  in  detection.  A  second  more  desirable  method 
would  be  to  rely  on  the  chi-square  test  for  detection.  Figure  F.9  shows  favorable  chi- 
square  results  in  terms  of  clean  and  quick  detection.  A  window  size  of  15  samples  will  give 
fast  detection  (two  seconds),  no  false  alarms  and  less  than  20  sec  of  delay  in  dropping  the 
failure  condition. 

The  next  problem  is  to  isolate  the  failure.  The  chi-square  test  provides  no  indication 
of  which  sensor  has  failed,  but  the  GLR  results  are  much  more  useful.  The  GLR  for  SVl 
has  a  very  large  value  at  2300  sec,  which  matches  the  time  of  the  dither  discussed  earlier, 
and  this  large  spike  distinguishes  it  from  the  other  matching  filters  (see  Figure  F.6).  A 
delay  of  almost  300  sec  would  exist,  which  is  unacceptable,  but  clear  identification  of  the 
failed  sensor  is  realized  based  on  a  second  threshold  designed  to  isolate  failures.  Further 
discussion  of  this  idea  will  be  presented  in  Chapter  V. 
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The  follow-on  goal  of  generating  corrective  feedback  to  the  Kalman  filter  was  not 
pursued  extensively,  but  a  quick  look  at  the  MLE  of  u  for  this  run  provided  some  interesting 
results.  Figure  F.8  shows  two  plots  of  v  with  different  vertical  scales.  The  actual  failure 
size  was  7000  ft,  and  close  inspection  of  these  plots  reveals  that  v  could  provide  a  crude 
estimate  of  the  size  of  the  failure  during  portions  of  the  failure  time.  The  same  problem  of 
poor  performance  prior  to  the  dither  signal  is  encountered  here,  along  with  disturbances  in 
v  after  3800  sec.  The  aircraft  dynamics  (an  altitude  change  at  3800  sec)  are  speculated  to 
have  caused  this  disturbance.  Continued  research  into  using  MLE’s  for  corrective  feedback 
is  encouraged  but  these  va'ues  must  be  checked  carefully  to  ensure  they  do  not  adversely 
affect,  the  filter  performance  by  providing  erroneous  feedback. 

Similar  results  were  obtained  for  run  five  in  which  the  bias  was  reduced  to  700  ft. 
A  single  GLR  plot  based  on  SV1  is  shown  in  Figure  F.10  along  with  its  corresponding 
fail  flag  result.  A  problem  was  encountered  with  missed  alarms  prior  to  the  dither  signal. 
Although  variations  in  the  threshold  value  and  number Jow  could  overcome  this  problem, 
these  variations  would  inhibit  detection  of  other  types  of  failures  discussed  previously.  The 
same  solution  to  use  the  chi-square  test  for  detection  was  available  and  it  was  determined 
that  a  threshold =50  would  be  effective  for  this  failure  (see  Figure  F.ll).  Other  threshold 
values  would  also  be  acceptable  for  this  failure,  but  using  a  threshold=50  coincided  with 
other  results.  Final  implementation  studies  could  be  accomplished  to  determine  the  best 
threshold  value  for  all  types  of  failures  being  considered. 

Ln  order  to  verify  the  influence  of  dither  signals  on  the  FD1  algorithms,  run  six  was 
conducted  with  the  failure  induced  during  the  highly  dynamic  time  frame  from  4000  - 
6000  sec.  Large  variations  in  the  GLR  and  x  data  correlated  to  quick  turns  in  the  flight 
profile,  as  anticipated.  The  GLR  plot  and  fail  flag  based  on  SVl  are  shown  in  Figures  F.1‘2 
and  F.13  with  unfavorable  results  occurring  with  nutnberJow—i\yo.  In  order  to  prevent 
missed  alaxms,  a,  nnmberJow=  45  was  required.  Again,  a  delay  occurred  in  dropping  the 
failure  condition,  but  the  chi-square  test  provided  the  necessary  alarm  with  the  same 
threshold =50  used  ea  lier.  The  same  technique  of  detection  through  the  chi-square  test 
and  isolation  via  the  GLR  test  would  work  1  his  failure. 
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The  last  bias  type  failure,  run  eight,  assumed  that  all  the  SV’s  would  be  spoofed 
simultaneously.  The  state  plots  in  Figures  F.15  -  F.17  show  the  robustness  of  the  filter  to 
adapt  to  this  failure.  Fewer  than  ten  updates  at  a  two-second  sample  rate  were  required 
for  the  filter  t  >  regain  accurate  tracking  of  all  states.  This  coincides  with  the  speculation 
in  Section  3.4.1  that  the  filter  would  attribute  this  failure  to  an  increase  in  the  user  clock 
bias.  Figure  F.18  shows  that  the  failure  is  quickly  detected  but  the  adaptive  nature  of  the 
EKF  prevents  a  sustained  recognition  of  the  failure.  Further  simulations  were  conducted 
with  various  biases  and  values  as  low  as  50  ft  were  easily  detected  by  the  GLR  algorithm. 
This  failure  was  not  pursued  further,  for  reasons  explained  in  Section  3.4.1. 

4-3.3  Ramp  Failures.  Runs  eight  and  nine  refer  to  ramp  tyne  failures  of  slopes  two 
and  one  ft/sec,  respectively.  The  state  plots  in  Figures  G.l  and  G.2  show  the  gradually 
increasing  degradation  in  the  filter  performance  due  to  this  type  of  spoofing.  The  residuals 
in  Figures  G.3  and  G.4  indicate  that  failure  detection  should  be  possible,  and  the  spikes 
in  Figure  G.3  highlight  the  effects  of  dithering  caused  by  high  dynamics. 

The  GLR  algorithm  was  somewhat  successful  in  detecting  the  larger  ramp  failure=2T 
with  some  delay.  Figure  G.5  illustrates  that  the  all  five  matching  filters  will  cause  missed 
and  false  alarms.  A  closer  look  at  these  results  is  shown  in  Figure  G.6  in  which  the  vertical 
scale  ha,s  been  reduced.  From  this  plot  the  matching  fillers  do  not  appear  to  provide  quick 
or  reliable  detection.  The  results  of  the  fail  flag  routine  are  shown  in  Figure  G.7  with  the 
thresholds  20  as  discussed  on  page  4-4  and  an  increase  in  number Jow~l0.  This  increase 
in  number-low  did  not  adversely  effect  previous  results  on  other  types  of  failures  other  than 
the  expected  delay  in  dropping  the  fail  flag.  The  two  major  problems  with  these  results 
are  the  356-sec  detection  delay  seen  in  even  the  best  matching  filter,  SV1,  and  the  multiple 
missed  and  false  alarms  seen  in  the  other  matching  filters.  A  positive  result  was  the  large 
spike  seen  at  the  time  of  the  dither  (2300  sec)  that  would  provide  isolation  of  tire  failure 
on  SV1. 

Again,  the  chi-square  test  provided  better  results  for  detection  with  the  ranip=-;2T. 
Figure  G.8  shows  reliable  detection,  based  on  a  threshold— 50,  beginning  250  sec  after  the 
failure  which  occurred  at  1500  sec.  The  250 -sec  delay  seems  excessive,  but  the  failure 
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is  subtle  and  the  aircraft  is  spoofed  approximately  400  ft  at  the  time  of  detection.  The 
chi-square  test  also  dropped  the  failure  more  quickly  than  the  GLR  test  and  did  not  suffer 
from  the  high  dynamics  between  5000  and  6000  sec. 

Results  with  the  ramp=lT  are  shown  in  Figures  G.9  -  G.ll.  The  GLR  test  was 
more  reliable  for  detection,  as  shown  in  Figure  G.10,  when  compared  to  the  ramp=2T 
results,  but  this  should  be  attributed  to  the.  lack  of  dynamics  during  the  failure  time 
frame.  The  chi-square  test  was  not  as  successful  with  this  more  subtle  failure  and  suffered 
the  same  delay  in  detecting  the  failure  as  did  the  GLR.  The  dither  signal  is  again  the  key 
to  detection,  and  isolation  was  also  possible  based  on  the  matching  fdter  results.  For  this 
particular  simulation  with  the  dither  signal  occurring  300  sec  after  the  failure,  the  aircraft 
would  only  be  spoofed  approximately  150  ft  off  course  prior  to  detection.  Although  these 
errors  due  to  detection  delay  might  be  acceptable  for  most  missions,  it  points  out  the  need 
to  conduct  a  study  on  optimal  input  signals  for  the  purpose  of  enhancing  failure  detection 
and  isolation.  Mehra  discusses  a  variety  of  techniques  for  determining  these  optimal  inputs 
including  observation  of  the  systems  eigenvalues  (18,  19).  Some  quick  and  simple  methods 
for  system  identification  are  presented  by  Zarrop  and  Goodwin  based  on  minimizing  a 
scalar  function  of  the  information  matrix  (32).  Consideration  should  also  be  given  to 
developing  failure  models  for  the  matching  filters  designed  to  detect  ramp  failures  with,  the 
expectation  that  detection  time  will  decrease. 

4-3.4  Failed  GPS.  A  complete  loss  of  one  pseudorange  entering  the  NllS  fdter  was 
induced,  and  the  degraded  performance  of  the  filter  is  shown  in  Figure  H.l.  The  filter  loses 
track  until  3800  sec  after  the  failure  is  removed,  revealing  a  major  liability  in  the  design 
of  the  NRS  filter.  Figures  II  .2  and  H.3  show  the  residuals  as  excessive  and  divergent, 
indicating  that  the  filter  is  unable  to  provide  a  good  navigation  solution  in  the  absence 
of  even  one  pseudorange  (PR).  One  solution  to  this  problem  is,  after  detecting  the  loss 
of  the  FR  causing  the  filter  divergence,  move  the  filter  back  in  time  by  recalculating  the 
Kalman  filter  equations  without  Lie  PR  measurement  input.  This  would  prevent  the  other 
measurements  (transponders,  baro-altimeter,  and  Doppler)  from  being  affected  by  the 
coupling  discussed  earlier,  and  the  filter  would  remain  stable.  Once  the  PR  input  was 
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reacquired,  this  measurement  would  be  processed  by  the  Kalman  filter  as  before.  This 
idea  was  not  pursued  due  to  time  constraints  but  is  mentioned  in  the  recommendations  for 
future  work.  The  FDI  results  are  skewed  by  the  filter  divergence  and  the  GLR  algorithm 
goes  unstable.  Defection  of  this  failure  is  not  a  major  concern  since  standard  GPS  receiver 
design  would  recognize  the  problem  and  compensate  if  possible.  Finally,  this  type  of  failure 
would  not  be  a  good  tactic  for  a  spoofer,  as  discussed  in  Section  3.4.1,  but  it  is  comforting 
to  know  that  this  failure  is  distinguishable  from  typical  spoofer -induced  failures. 

4-3.5  Numerical  Precision  and  Modelling.  Considerable  attention  was  given  to  the 
numerical  precision  of  the  computers  running  the  simulations  and  the  FDI  algorithms.  In¬ 
vestigation  into  the  matrices  passed  from  MSOFE  to  Matrixx  revealed  a  potential  problem. 
The  two  rows  of  the  measurement  matrix,  H,  corresponding  to  states  one  and  two  were  10s 
larger  than  many  of  the  other  rows,  but  both  computing  systems  were  using  double  preci¬ 
sion  in  all  their  calculations,  so  there  might  not  be  a  problem.  The  positive  performance 
of  both  the  Kalman  filter  within  MSOFE  and  the  GLR  algorithm  in  Mat,rixx  indicated 
that  everything  was  fine,  but  it  is  possible  that  their  performance  could  be  improved  by 
performing  a  similarity  transformation  on  the  filter  model  to  scale  the  states  more  evenly. 
Details  for  this  type  of  transformation  are  showm  by  Maybeck  (13:p.  28).  Time  limitations 
prevented  this  work  from  being  completed. 

Another  factor  affecting  the  performance  of  the  GLR  test  is  mismodelling  even  though 
the  NRS  filter  was  shown  to  exceed  specifications  in  providing  a.  navigation  solution.  The 
dependence  of  the  GLR  algorithm  on  the  system  dynamics  makes  it  susceptible  to  inaccu¬ 
racies  in  modelling,  especially  when  the  truth  model  has  97  states  versus  the  filter  model  of 
only  15  states.  A  quick  study  was  performed  in  which  both  the  truth  and  filter  models  had 
only  15  states  to  ensure  that  modelling  was  not  a  major  factor  in  the  GLR  performance. 
No  noticeable  improvement  resulted  from  this  study.  Detection  delays  on  subtle  failures 
were  not  significantly  reduced  and  failure  isolation  was  still  not  possible  prior  to  the  dither 
signals.  These  findings  indicated  that  higher  order  models  for  the  filter  would  not  improve 
FDI  performance  sufficiently  to  justify  the  increase  in  computational  loading  required  for 
a  larger  filter. 
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4-4  Summary 

This  chapter  focussed  on  the  results  obtained  in  the  research.  A  brief  discussion 
of  the  piecewise-constant  nature  of  the  dynamics  matrix  and  concerns  about  numerical 
precision  and  modelling  indicated  that  these  assumptions  were  valid.  Results  on  filter 
performance  were  very  good  and  a  comparison  to  past  research  showed  that  the  reduced 
order  filter  model  worked  well.  Analysis  of  the  FD1  performance  revealed  some  positive 
results  but  techniques  for  enhancing  this  performance  would  be  necessary  prior  to  final 
implementation  on  actual  aircraft.  Chapter  V  offers  some  suggestions  to  accomplish  this 
goal  and  provides  some  ideas  for  future  research. 
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V.  Conclusions  and  Recommendations 


This  chapter  provides  a  brief  summary  of  the  results,  including  a  possible  method 
of  implementation  for  the  FDI  system.  Recommendations  for  future  research  are  also 
presented. 

5.1  NRS  Filter  Performance 

The  15-state  NRS  filter  exceeded  all  operational  specifications  in  its  ability  to  pro¬ 
vide  an  accurate  navigation  solution.  Although  the  filter’s  performance  was  degraded  as 
compared  to  higher  order  models  used  in  prior  research,  the  state  estimates  were  reliable 
and  accurate.  Increases  in  tracking  errors  were  experienced  following  a  series  of  harsh 
maneuvers  by  the  aircraft,  prompting  the  use  of  adaptive  tuning  techniques  on  states  one, 
two,  12  and  14.  A  simple  technique  required  these  states  to  have  lower  and  upper  process 
noise  values  for  adaptation. 

The  reduced  order  RRS  filter  model  provided  adequate  estimation  for  the  transpon¬ 
ders  while  significantly  decreasing  the  number  of  states  in  the  filter  from  26  to  two.  The 
importance  of  this  simplification  is  ihe  reduction  in  the  load  on  the  computer  tasked  to 
run  the  Kalman  filter  algorithm.  This  accomplishment  is  magnified  by  the  fact  that  many 
of  the  small  aircraft  on-board  computers  in  the  Air  Force  inventory  lack  the  capability  to 
run  a  Kalman  filter  with  over  60  states  and  a  two-second  update  rate. 

5.2  FDI  Performance 

The  discretization  process  used  on  the  state  transition  matrix  assumed  time  invari¬ 
ance  over  the  two-second  sample  period  and  was  sufficient  to  generate  accurate  inputs 
for  the  GLR  algorithm.  This  resulted  from  the  combination  of  the  two-second  sampling 
period  and  the  8th  order  Pade  approximation  for  the  matrix  exponential  function.  Longer 
sample  periods  or  cruder  approximation  techniques  might  cause  the  discretization  process 
to  be  degraded,  resulting  in  improper  modelling  within  the  GLR  algorithm. 

A  combination  of  the  GLR  and  chi-squa/re  algorithms  provided  adequate  detection 
and  isolation  for  most  of  the  failures  considered.  Heavy  and  medium  jamming  levels  were 
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reliably  detected  but  light  jamming  was  not  noticable  above  the  noise  floor.  Although 
jamming  detection  using  the  FDI  algorithms  was  desired,  many  other  forms  of  detection 
are  available  using  current  electronic  warfare  techniques.  Even  more  important  is  the 
distinction  between  the  outputs  of  the  GLR  matching  filters  when  compared  to  each  other 
and  also  when  compared  to  other  types  of  failures  such  as  spoofing.  It  is  necessary  for  the 
FDI  system  to  recognize  the  type  of  failure  prior  to  generating  any  corrective  feedback  into 
the  NRS  filter  and  to  alarm  the  aircraft  pilot  properly  of  the  hostile  environment. 

The  FDI  algorithms  performed  well  for  both  types  of  spoofing.  Thresholds  were 
found  empirically  with  substantial  tradeoffs  required  to  prevent  false  and  missed  alarms. 
All  bias  values  were  easily  detected  and  isolated  using  a  combination  of  the  chi-square 
test  for  detection  and  the  GLR  test  for  isolation.  The  chi-square  test  minimized  detection 
delay  to  two  sec,  eliminated  false  alarms  a,nd  minimized  the  delay  in  dropping  the  failure 
condition  to  20  sec.  If  the  FDI  algorithms  were  attempting  to  detect  failures  in  a  flight 
control  system  rather  than  a  navigation  system,  then  these  delays  would  be  unacceptable, 
but  for  the  scenario  described  in  Chapter  I,  these  delays  should  not  cause  major  problems 
in  completing  mission  objectives.  Ramp  failures  presented  additional  problems  for  the 
chi-square  test.  The  gradual  changes  in  the  residuals  caused  delays  in  detection  with 
larger  ramp  values  having  delays  of  about  four  minutes  and  small  ramp  values  completely 
avoiding  detection  prior  to  dithering  the  filter.  These  results  would  not  be  acceptable  for 
final  implementation. 

Two  important  aspects  must  be  addressed  for  implementation  of  a  reliable  FDI  sys¬ 
tem.  First,  a  two-level  scheme  should  be  used  in  which  the  chi-square  test  is  utilized  for 
detection  and  the  GLR  test  provides  isolation.  It  should  be  realized  that  both  algorithms 
use  the  inverse  of  the  residual  covariance,  in  their  calculations,  so  the  chi-square 

test  does  not  generate  much  computational  overhead  assuming  the  GLR  test  is  already 
being  conducted.  The  GLR  algorithm  does  require  a  significant  number  of  computations 
and  may  interfere  with  chi  -square  calculations  done  by  an  on  board  computer,  but  the 
two-level  design  allows  for  priority  to  be  given  to  the  chi-square  test  for  detection  while 
the  GLR  test  can  run  on  a  lower  priority  lor  isolation.  A  study  on  this  idea  could  provide 
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optimal  performance  of  the  overall  FDI  scheme.  Neither  algorithm  could  independently 
outperform  this  two-level  scheme  in  terms  of  detection  and  isolation. 

The  second  major  concern  for  implementation  is  whether  purposely  dithering  the 
system  could  enhance  FDI  performance.  Results  clearly  indicate  that  harsh  aircraft  ma¬ 
neuvers  enhanced  isolation  for  the  GLR  test  on  all  failures  and  detection  for  the  chi-square 
test  on  subtle  or  small  failures.  Maybeck  discusses  the  use  of  probing  with  the  intention 
that  it  “ purposely  excite  certain  modes  of  the  system  in  order  to  aid  the  identification  of 
uncertain  parameters”  (15:p.  229).  One  tactic  discussed  is  the  use  of  “S-turns”  performed 
by  the  pilot  and  a  simple  variation  in  the  flight  profile  would  reveal  the  usefulness  of  this 
idea.  Stratton,  Menke  and  Hanlon  looked  into  continuous  and  periodic  dithering  signals 
induced  through  control  inputs  (6,  20,  28).  This  would  change  Equation  (2.1)  to 

x(t)  =  f  [x(i),  <,u(<)]  +  G(t)w(t)  (5.1) 

where  u(t)  —  control  input,  A  study  could  be  conducted  on  optima)  control  inputs  for  this 
navigation  system  to  determine  the  best  choices  for  this  method.  Preliminary  indications 
are  that  rather  large  dither  signals  are  required  for  this  system  so  pilot-induced  dithers 
would  probably  work  better  than  automatically  induced  subliminal  probes  or  automatically 
induced  nonsubliminal  probes  that  would  bounce  around  the  aircraft  and  its  crew. 

Comparison  of  a  15-state  filter  model  against  a  15-state  truth  model  showed  that 
higher  order  filter  models  would  not  significantly  enhance  FDI  performance.  Numerical 
precision  was  not  considered  to  cause  a  major  problem  in  the  performance  of  the  Kalman 
filter  or  the  FDI  algorithm,  although  repealing  the  states  should  still  be  considered. 

5.2.1  Corrective  Feedback.  Although  time  constraints  prevented  research  into  feed¬ 
back  techniques  designed  to  correct  for  failures,  some  observations  can  be  made.  As  men¬ 
tioned  earlier,  adaptive  tuning  techniques  in  which  Q  was  increased  were  employed  to 
ensure  adequate  filter  performance  under  normal  conditions.  This  concept  could  be  ex¬ 
tended  to  enhance  filter  performance  once  a  failure  was  detected.  The  FDI  system  is  i  aw 
confronted  with  a  third  task  beyond  that  of  detection  and  isolation.  Assuming  a  bias  failure 
has  been  induced,  the  FDI  system  must  estimate  the  size  of  the  bias  in  order  to  determine 


the  amount  of  adaptive  feedback  (represented  by  changes  in  Q)  needed  to  compensate  for 
this  failure.  If  the  estimate  is  too  large,  the  filter  may  become  too  conservatively  tuned  and 
performance  will  degrade.  If  the  estimate  is  too  small,  insufficient  tuning  may  occur  and 
the  filter  will  not  be  able  to  track  the  states  closely.  Additional  problems  may  result  from 
the  filter  becoming  accustomed  to  the  errors  induced  by  the  failure.  If  the  filter  adapts 
to  the  bias,  it  may  become  less  aware  of  minor  changes  in  the  bias  and  might  lose  the 
ability  to  determine  when  the  failure  has  been  removed.  Since  the  pilot  would  be  alarmed 
of  the  initial  failure,  operator-induced  dither  signals  could  be  used  to  help  improve  future 
detection  and  the  overall  filter  performance  would  be  improved 

Another  method  of  compensation  is  to  estimate  the  size  of  the  failure  and  remove  it 
from  the  incoming  measurement.  Assume  the  failure  appears  as  a  bias,  b,  on  the  measure¬ 
ment  equation  taken  from  Equation  (2.11). 

Szn{ti)  =  K[ti-,xn{ti)]6xn(t)  +  v(i<)  +  b  (5.2) 

If  the  GLR  algorithm  were  able  to  find  0  accurately  which  is  the  MLE  of  the  size  of  the 
failure,  then  simply  subtracting  0  from  Equation  (5.2)  would  eliminate  the  failure  in  the 
system  and  the  output  of  the  Kalman  filter  would  not  be  degraded.  Other  algorithms  such 
as  MMAE  could  also  use  this  idea. 

Recommendations 

A  brief  list  of  recommendations  for  future  research  is  presented,  with  many  of  the 
details  concerning  these  items  presented  earlier. 

•  1  inaiize  the  similarity  transform  to  rescale  the  states  and  compare  results  to  those 
achieved  without  rescaling  for  both  the  Kalman  filter  and  FD1  algorithms. 

•  Verify  dominance  of  the  diagonal  terms  in  the  residual  covariance  matrix,  A  (  /.;),  over 
the  non -diagonal  terms  of  this  matrix.  See  assumption  7. 
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Perform  a  study  on  optimal  inputs  intended  to  dither  the  system.  Improvements 
should  be  noted  in  the  NRS  filter  performance,  but  the  main  goal  would  be  to 
decrease  delays  in  detection  and  isolation  of  failures. 

Pursue  techniques  for  feeding  back  corrective  signals  such  as  adaptive  tuning  or 
measurement  correction  discussed  previously.  Careful  investigation  of  errors  in  the 
state  estimates  should  indicate  the  success  of  these  methods. 


Consider  changing  the  failure  models  for  the  matching  filters  to  look  for  ramp  failures 
rather  than  biases.  A  scheme  with  10  matching  filters  (five  for  biases  AND  five  for 
ramps)  should  be  compared  against  those  with  only  five  matching  filters  (five  for 
biases  OR  five  for  ramps). 

Look  into  theoretical  techniques  for  establishing  thresholds  rather  than  finding  them 
empirically. 

Change  the  GLR  algorithm  listed  in  Appendix  I  to  allow  for  various  failure  times. 
This  implies  that  0  could  take  on  any  value  within  the  window,  resulting  in  an  MLE 
of  6  and  a  GLR  as  a  function  of  t{  and  §.  This  GLR  should  better  identify  the  time  of 
the  failure  and  may  outperform  the  simplified  GLR  algorithm.  Consideration  would 
have  to  be  given  to  increased  detection  time  caused  by  the  increased  computations 
necessary  to  implement  this  new  GLR  test. 


Look  into  MMAE  techniques  as  a  replacement  for  the  EDI  scheme.  MMAE  could 
oe  used,  m  conjunction  wim  many  ui  me  recoiiuueiidatiGiis  listed  above  and  the  chi 
square  test  should  always  be  considered  as  an  additional  source  of  information.  The 
most  significant  task  associated  with  this  idea  is  developing  the  software  capable  of 
performing  the  simulations.  MSOFE  is  not  currently  structured  to  handle  the  bank 
of  Kalman  filters  associated  with  an  MMAE  scheme.  The  alternative  would  be  to 
adopt  another  software  package  designed  to  run  multiple  Kalman  filters.  The  major 
problem  would  be  to  convert  the  models  and  particularly  the  overhead  associated  with 
generating  the  measurements  for  the  GPS  and  RRS  portions  of  the  NRS  system.  In 
either  case,  a  head-to-head  comparison  of  the  results  achieved  by  the  MMAE  method 


and  those  achieved  in  this  thesis  should  provide  insights  into  the  relative  strengths 
and  weaknesses  of  the  GLR  and  MMAE  approaches. 
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Appendix  A.  Error  M:del  State  Definitions 


Tabular  listings  of  the  truth  and  filter  models  are  presented.  Tables  A.l  and  A.2 
show  the  41-state  INS  truth  model  with  the  LN-93  state  numbers  given  for  reference  to 
the  Litton  technical  report  on  the  INS  (9).  Tables  A. 3  and  A.4  list  the  RR.S  and  GPS 
states  respectively  and  Table  A. 5  lists  only  the  states  used  in  the  NBS  filter  model. 
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Table  A.l.  41-State  INS  System  Model:  First  20  States 


State 

Number 

State 

Symbol 

Definition 

H29 

1 

69, 

X-component  of  vector  angle  from  true  to  computer  frame 

i 

2 

kiA 

Y-component  of  vector  angle  from  true  to  computer  frame 

2 

3 

69, 

Z-component  of  vector  angle  from  true  to  computer  frame 

3 

4 

<t>x 

X-component  of  vector  angle  from  true  to  platform  frame 

4 

j  5 

(j>y 

Y-component  of  vector  angle  from  true  to  platform  frame 

5 

j  6 

<t>z 

Z-component  of  vector  angle  from  true  to  platform  frame 

6 

7 

svx 

X-component  of  error  in  computed  velocity 

7 

8 

msm 

Y-component  of  error  in  computed  velocity 

8 

9 

6V, 

Z-component  of  error  in  computed  velocity 

9 

10 

6h 

Error  in  vehicle  altitude  above  reference  ellipsoid 

10 

11 

bhB 

Total  baro-altimeter  correlated  error 

23 

16 

ShL 

Error  in  lagged  inertial  altitude 

SH 

17 

ss3 

Error  in  vertical  channel  aiding  state 

SB 

mm 

■S 

Error  in  vertical  channel  aiding  state 

13 

19 

m 

X-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

17 

20 

m 

Y-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

18 

21 

*>c 

Z-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

19 

22 

tgx 

X-component  of  gravity  vector  errors 

20 

23 

% 

Y-component  of  gravity  vector  errors 

21 

24 

6g* 

Z-component  of  gravity  vector  errors 
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Table  A. 2.  41- State  INS  System  Model:  Second  21  States 


State 

Number 

Definition 

l,N-93 

State 

25 

bx 

X-component  of  gyro  drift  rate  repeatability 

30 

26 

K 

Y-component  of  gyro  drift  rate  repeatability 

31 

27 

bz 

Z-component  of  gyro  drift  rate  repeatability 

32 

28 

mm 

X-component  of  gyro  scale  factor  error 

33 

j  29 

mm 

Y-component  of  gyro  scale  factor  error 

34 

30 

Ss, 

Z-component  of  gyro  scale  factor  error 

35 

■SB 

vbt 

X-component  of  accelerometer  bias  repeatability 

48 

SB 

1 

Y-component  of  accelerometer  bias  repeatability 

49 

i  -i** 

Ul# 

mm 

Z-component  of  accelerometer  bias  repeatability 

50 

34 

m 

X-component  of  accelerometer  and  velocity- 
quantizer  scale  factor  error 

51 

35 

m 

Y-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 

52 

36 

m 

Z-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 

53 

37 

Sqa* 

X-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 

54 

38 

Sqav 

Y-component  of  accelerometer  and  velocity 

- j.: - - 1  _  _ 

qudilitiZitll  labiui  aoj  imucmjf 

55 

39 

Sqa. 

Z-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 

56 

40 

Mi 

X  accelerometer  misalignment  about  Z-axis 

66 

41 

M2 

Y  accelerometer  misalignment  about  Z-axis 

67 

42 

Ms 

Z  accelerometer  misalignment  about  Y-axis 

68 

43 

<7l 

X-accelerometer  misalignment  about  Y-axis 

69 

44 

<72 

Y-accelerometer  misalignment  about  X-axis 

69 

45 

<73 

Z-accelerometer  misalignment  about  X-axis 

69 
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Table  A.3.  20-State  RItS  System  Model 


State 

Number 

| 

Definition 

12 

Sftt 

Range  error  due  to  equipment  bias 

13 

Svb 

Velocity  error  due  to  equipment  bias 

46 

M'ti* 

Transponder  1  x-component  of  position  error 

47 

SPi  i. 

Transponder  1  y-component  of  position  error 

48 

fPT\. 

Transponder  1  z-component  of  position  error 

mm 

6RTl„ 

Transponder  1  range  error  due  to  atm  propagation 

SB 

mm 

Transponder  2  x-component  of  position  error 

51 

o 

c: 

Transponder  2  y-component  of  position  error 

52 

&PT2, 

Transponder  2  u-componcnt  of  position  error 

53 

fjP.T2r 

54 

$Pr*x 

Transponder  3  x-component  of  position  error 

mmm 

t>PT3„ 

Transponder  3  y-component  of  position  error 

56 

6Prs. 

Transponder  3  z-component  of  position  error 

57 

6  Hx3a 

Transponder  3  range  error  due  to  atm  propagation 

58 

6Pr*. 

Transponder  4  x  component  cf  position  error 

59 

SPta, 

Transponder  4  y-component  of  position  error 

60 

SPta, 

Transponder  4  z-component  of  position  error 

61 

SPj‘4. 

62 

6Pn. 

Transponder  5  x-component  of  position  error 

62 

8h^. 

Transponder  5  y -component  of  position  error 

6J 

6Prs. 

Transponder  5  z-component  of  position  error 

65 

ftPrb, 

Transponder  5  range  error  due  to  atm  propagation 

66 

fiPrc* 

Transponder  6  x- component  of  position  c;ror 

67 

6Pts. 

Transponder  6  y-componcnl  of  position  error 

68 

SPt*. 

Transponder  6  z-component  of  position  error 

69 

SR-T  6„ 

Transponder  6  range  error  due  to  atm  propagation 

AM 


Table  A.4.  30-State  GPS  System  Model 


Defini  tion 


User  clock  bias 
User  clock  drift 
SV  1  code  loop  error 
SV  1  tropospheric  error 
SV  1  ionospheric  error 
SV  1  clock  error 

SV  1  x-component  of  position  error 
SV  1  y-component  of  position  error 
SV  1  zy component  of  position  error 
SV  2  code  loop  error 
SV  2  tropospheric  error 
SV  2  ionospheric  error 
SV  2  clock  error 

SV  2  x-component  of  position  error 
SV  2  y-component  of  position  error 
SV  2  z-component  of  position  error 
SV  3  code  loop  error 
SV  3  tropospheric,  error 
SV  3  ionospheric  error 
SV  3  clock  error 

SV  3  x-component  of  position  error 
SV  3  y-component  of  position  error 
SV  3  z-component  of  position  error 
SV  4  code  loop  error 
SV  4  tropospheric  error 
SV  4  ionospheric  error 
SV  4  clock  error 

SV  4  x-componcnt  of  position  error 
SV  4  y-component  of  position  error 
SV  4  z-component  of  position  error 


Table  A. 5.  15-State  Reduced-Order  Filter  Model 


A-G 


Appendix  B.  Dynamics  Matrices  and  Noise  Values 


B.l  Definition  of  Dynamics  Matrices 

The  LN-93  error-state  dynamics  matrix  F  is  defined  in  Chapter  III  as  a  combination 
of  submatrices.  The  NON-ZERO  elements  of  these  submatrices  are  presented  in  the  tables 
which  follow.  All  the  variables  shown  in  the  following  tables  are  defined  in  the  LN-93 
technical  report  along  with  their  unns  (9). 


Table  B.l.  Elements  of  the  Dynamics  Submatrix  Fu 


Table  B.3.  Elements  of  the  Dynamics  Submatrix  Fj3 


Term 

Element 

Term 

Element  Tei  m 

(4,25) 

warn 

(4,26) 

C12 

(4,27) 

r  r 

^13 

(4,28) 

(4,29) 

I:  Mil 

(4,30) 

(5,25) 

Ai 

(5,26) 

C22 

(5,27) 

E23 

(5,28) 

(5,29) 

(5,30) 

C23  Wii, 

(6,25) 

C31 

1  ■£&$■ 

r 

'-'32 

(6,27) 

— 

(6,28) 

(6,29) 

Bi 

(6,30) 

Table  B.4.  Elements  of  the  Dynamics  Submatrix  Fm 


|  Element 

Term 

— 

Element  J  Term  ] 

Element 

Term 

(7,31) 

Ai 

Emm 

(7,33) 

£*13 

(7,34) 

warnm 

(7,35) 

Esmm 

(7,36) 

CuA? 

(7,37) 

Cn\A°\ 

wsm m 

C\M*\ 

(7,39) 

CjzjA?’  1 

A3  A® 

(T ,40) 

(7,41) 

mam 

(7,42) 

(7,43) 

(8,31) 

C21 

(8,32) 

A2 

(8,33) 

<^23 

(8,34) 

AiA? 

(8,35) 

A2A? 

(8,36) 

As^f 

(8,37) 

Ai|ab| 

(8,38) 

A2|A?| 

(8,33) 

C*\A?'\ 

(8,40) 

(8,41) 

-AaA? 

(8,42) 

(8,43) 

(9,31) 

Ai 

(9,32) 

A  2 

(9,33) 

Qb 

(9,34) 

Emm 

(9,36) 

AsAf' 

(9,37) 

Ai|A®|  1 

(9,38) 

A2I  Af  1 

(9,39) 

A3M?  I 

(9,40) 

SUP 

(9,41) 

(9,42) 

A.V^' 

(9,43) 

Bssi 

li-3 


Table  B.5.  Elements  of  the  Dynamics  Submatrix  F2? 


Element 

Term 

Element 

Term 

Element 

Term 

(19,19) 

191 

(21,21) 

"fa-* 

(22,22) 

mm 

(24,24) 

(11,11) 

mm 

II 

il^CTi 

The  NON-ZERO  elements  of  the  dynamics  matrix  representing  the  the  GPS  and 
RRS  are  shown  in  the  following  table. 


Table  B.8.  Elements  of  the  Dynamics  Matrix  for  GPS  &  RRS 


(49,49) 

-1/300  ft2 /sec 

(53,53) 

-1/300  ft2 /sec 

(57,57) 

-1/300  ft2 /sec 

(61,61) 

-1/300  ft2 /sec 

(65,65) 

-1/300  ft2 /sec 

(69,69) 

-1/300  ft2 /sec 

(70,70) 

-1  ft2 /sec 

(71,71) 

-1/500  ft2 /sec 

(72,72) 

-1/1500  ft2 /sec 

(77,77) 

-1  ft2/  sec 

(78,78) 

-1/500  ft2 /sec 

(79,79) 

-1/1500  ft2 /see 

(84,84) 

-1  ft2  /sec 

(85,85) 

-1/500  ft2 /sec 

(86,86) 

-1/1500  ft2 /sec. 

(91,91) 

-1  ft2 /sec 

(92,92) 

-1/500  ft2 /sec 

(93,93) 

-1/1500  ft2  /sec 

!  (14,16) 

1  ft2 /see 

B.2  Elements  of  the  Process  Noise  and  Measurement  Noise  Matrices 

The  process  noise  strength  matrix  Q  associated  with  the  INS  truth  model  is  also 
partitioned  into  submatrices  as  described  in  Chapter  111.  The  NON-ZERO  elements  of 
these  submatrices  are  shown  in  Tables  B.7  and  B.8.  Note  that  the  cr2  term*  m  these 
two  tables  are  variable  names  only  as  defined  in  the  Litton  technical  report  2nd  do  not 
represent  variance  terms  typically  associated  with  cr2.  The  process  noise  for  the  GPS  and 
HRS  portions  of  the  truth  model  are  listed  in  Table  U.9.  Finally,  the  process  noise  values 
used  in  the  filter  and  the  measurement  noise  values  R  are  presented  in  Tables  B.10  ana 
B.ll. 


li-4 


Table  B.7.  Elements  of  Truth  Model  Process  Noise  Submatrix  Qn 


Element 

Term 

Element 

Term 

(4,4) 

<r 

(5,5) 

<„ 

(6,6) 

< 

(7J) 

< 

(8,8) 

< 

(9,9) 

Table  B.8.  Elements  of  Truth  Model  Process  Noise  Submatrix  Q22 


|  Element 

Term 

Element 

Term 

Element 

Term 

1  (H,ll) 

2/?vVc<4„. 

msm 

2ftsxCTL 

MEMIB 

II  ■ 

Table  B.9,  Elements  of  Truth  Model  Process  Noise  for  GPS  &.  RRS 


(49,49) 

6.667x1  O'13  ft2 /sec 

I1I2&S2I 

6.667xl0-13  ft2 /sec 

- 1 

Cn 

J-4 

"cn 

-4 

6.667xl0-13  ft2 /sec 

6.667xl0~13  ft3 /sec 

(65,65) 

6.667xl0-13  ft2 /sec 

(69,69) 

6.667xl0~13  ft2 /sec 

(  7(\  7^ 

V  .  V,  . 

0.5  ft2 /sec 

(71,71) 

0,004  ft3 /sec 

(72,72) 

0.004  ft2 /sec 

(77,77) 

0.5  ft2 /sec 

(78,78) 

0.004  ft2 /sec 

(79,79) 

(84,84) 

0.5  ft'2 /sec 

(85,85) 

(86,86) 

0.004  ft2 /sec 

(91,91) 

0.5  ft2 /sec 

(92,92) 

0.004  ft2 /sec 

(93,93) 

Table  B.10.  Filter  Process  Noise  Q 


Term 

Element 

Term 

■Ha 

O.lxlO-13  &  l.OxlO-13  rad2 /sec 

(2,2) 

0.58xl0-13  &  l.OxlO-13  rad2 /sec 

(3,3) 

0.0  rad2 /sec 

■CBM 

100  rad2 /sec 

(5,5) 

500  rad2 /sec 

KaaH 

45  rad2 / sec 

(7,7) 

800  ft2 /sec3 

(8,8) 

(9,9) 

(10,10) 

(11,11) 

400  ft2 /sec2 

(12,12) 

17  &  40  ft2 /sec2 

[i  (13,13) 

0.0  ft2 /sec3 

(14,14) 

14  &  40  ft2 /sec2 

(15,15) 

0.5x10" 13  ft2 /see? 

■■■ 

Table  B.li.  Truth  and  Filter  Measurement  Noises  R 


|  ideasurement 

Truth  Noise  j  Filter  Noise 

||  Baro  Altimeter 

r  2500  ft 2  J  2500  ft 2 

Appendix  C.  Miscellaneous  Plots 
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(degrees)  Latitude  (degrees) 


Figure  C.l.  Two  Hour  Flight  Profile 


C-2 


Figure  C.3.  Spoofer  Induced  Errors 


Appendix  D.  Baseline  Filter  Plots 


All  the  state  plots  contained  in  this  appendix  contain  five  traces.  The  innermost 

trace  ( - )  on  each  data  plot  is  the  mean  error  time  history  for  the  applicable  state. 

Mean  Error  is  defined  as  being  the  difference  between  the  filter’s  estimate  of  the  state  and 
the  true  state,  averaged  over  the  number  of  Monte  Carlo  runs  performed.  The  equation 
describing  this  relationship  is  defined  by  (13,  27): 


Md(ti)  —  i)  —  Jq  Y  Xiruej(ti)}  (^d) 

V  j- i  J=i 

where  Xj(U)  is  the  filter-computed  estimate  of  a  given  state  and  xtruej(ti)  is  the  truth 
model  value  of  the  same  state,  at  time  for  run  j,  and  N  is  the  number  of  time  histories 
in  the  simulation  (10  in  this  thesis). 

In  addition  to  the  center  trace,  two  more  puirs  of  traces  are  plotted  and  labeled 
Mean+-Sigrna.  The  first  pair  (represented  by  •••)  is  symmetrically  displaced  about  the 
mean  and  as  a  result  follows  the  “undulations”  of  the  Mc{t o-  ,|  !'e  locus  of  these  traces  is 
calculated  from  Me(U)  ±  >  where  PC(U)  is  the  true  ei.  variance  at  time  t,.  The 

true  standard  deviation  is  calculated  from  (13, 


Y  '  -J—MB.U)  \D.2) 

r-  i  v  —  i 

where  N  is  the  number  of  runs  in  the  Monte  Carlo  simulation  (10  in  this  thesis),  and  M,2(<<) 
is  the  square  of  the  mean  of  a  given  state  at  each  time  interest  (such  as  measurement 
times). 

The  last  pair  of  traces  ( - )  represent  the  filter  computed  ±  crjiHeT  values  for  the 

same  states  and  are  symmetrically  displaced  about  zero  because  the  filter  “believes”  that 
it  is  producing  zero-mean  errors  (15,  27).  These  quantities  are  propagated  and  updated  in 
the  MSOFE  (22,  27)  software  using  the  covariance  propagation  equation  shown  in  Chapter 
II.  These  traces  represent  the  filter’s  estimate  of  its  own  error. 


O  fr 


»e(0  =  = 


IV 
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Similar  statistics  are  computed  for  the  measurement  residuals  with  the  residual  co- 
variance  defined  in  Chapter  II  and  labeled  Filter+Sigrua.  Finally,  baseline  plots  for  the 
GLR.  and  chi-square  algorithms  are  shown,  with  three  window  sizes  presented  for  the 
chi-square  test. 
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<1)  Latitude  Error  (Ft) 
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Figure  D.l. 
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Figure  D.3.  State  Plots;  Baseline 
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(4)  E  Tilt  Error  (arcEec>  (4>  E  Tilt  Error  (areBec) 


Figure  D.4.  Stale  Plots:  Baseline 
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(5)  N  Tilt,  Error  (arcaec)  (S)  N  Tilt  Error  (ftrcBec) 
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Figure  D.5.  State  Plots:  Baseline 
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(6)  Azimuth  Error  (arcsec)  (6)  Azimuth  Err*o  r 
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Figure  D.6.  State  Plots:  Baseline 
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Figure  1).7.  State  Plots:  Baseline 
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Plots:  Baseline 
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(12)  RRS  Range  Bias  ( ft ) 


RRS  Vel  Bias  (fpo)  (13)  RRS  Ve)  Bias  (fps) 


.06 


Mean  error 
Meanf-Sigma 


.04 


.02 

0 

-.02 


-  Filtert-Sigma 

■l  1  i  i  i  1  i  i  l  i.l  i,  ,-i  ,  i  Juj  i — i  i  i  1  i  i  1,1,1 — L— i,  .  i  — :  —i 

1000  2000  3000  4000  5000  6000  7000 

Time  (see) 


8000 


(b) 


Fi 


0  20 


CIS)  User  Clock  Drift  Err  (ft 


AS 


sv 


30 
20 
10 
0 

-10 

-20 

-30 

0  1000  2000  3000  4000  5000  6000  7000  6000 


Residual 

Residual+-Sigma 


l — * 


-  Filter+^kjma 

1  1  '  1  I  '  '  '  1  I  '  1  '  1  I  1  ' _ 1 _ '  I  ‘  1  '  ' —  I  .  ■  ■■ _ ] _ 1  1  I  „L.  | _ i  i  i _ L 


Time  (sec) 


* 

> 

Cfl 


20 
10 
0 

-10 
-20 
-2D 

V  0  1000  2000  3000  4000  5000  6000  7000  6000 


Residual 

Residual+-Sigma 


-  Filtert-Sigma 

J_l - 1 _ I  I  I,  _) - 1  ‘  '  '  1  I  '  '  '  '  I _ ...  .  _ I _ I _ I _ I _ I _ I _ I  -l  l  ■  -L  1  '  ■  .  I 


Time  (sec) 


Figure  D.17.  Residual  Plots:  Baseline 
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Figure  D.18.  Residual  Plots:  Baseline 
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Figure  D.19.  Residual  Plots:  Baseline 
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Appendix  F.  Jamming  Failure  Plots 


Selected  state  and  residual  plots  are  shown  to  illustrate  the  performance  of  the  NRS 
filter  under  heavy  jamming  conditions.  GLR  and  CIII  plots  are  presented  for  all  three 
jamming  levels. 
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Figure  E.l.  State  Plots:  Heavy  Jamming 
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(12)  RRS  Range  Bias  (ft)  (lO)  Altitude  Error-  ( ft) 


Figure  E.2.  State  Plots:  Heavy  Jamming 
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Figure  E.7.  GLlt:  Heavy  Jamming 
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Figure  E.8,  Fail  Flag:  Heavy  Jamming 
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Appendix  F.  Bias  Failure  Plots 


Selected  state  and  residual  plots  are  shown  to  illustrate  the  performance  of  the  NRS 
filter  with  bias  failures  induced  on  the  SV’s.  Some  of  the  statistics  are  masked  by  the  large 
vertical  scale  cf  the  plots,  but  this  allows  the  full  range  of  the  mean  error  and  residual 
values  to  be  seen,  which  best  characterizes  the  filter  performance.  A  variety  of  GLR  and 
X  plots  illustrate  the  FDI  performance  for  all  bias  failures  considered.  Several  plots  are 
shown  twice  with  two  different  vertical  scales  to  reveal  significant  details  that  are  lost  in 
other  plots. 
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Figure  F.5.  GLR:  Run  4,  Bias=7000 
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Figure  F.7.  Fail  Flag:  Hun  4,  Bias=7000 
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Figure  F.8.  MLE  £(4):  Run  4,  Bias=7000 
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Figure  F.9.  CHI:  Run  4,  Bias=7000 
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Figure  F.FJ.  GLR:  Run  6,  Buis -700,  4000-6000  sec 
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Figure  F.13.  Fail  Flag:  Run  6,  Bias=700,  4000-6000  see 
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Figure  F.14.  CHI:  Run  G,  Bias=700,  4000-6000  sec 
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Figure  F.1,5.  State  Plots:  Run  7,  Bias=700  on  All 
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Figure  F.18.  GLR/Fail  Flag:  llun  7,  Bias=700  on  All 
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Appendix  G.  Ramp  Failure  Plots 


Selected  state  and  residual  plots  are  shown  to  illustrate  the  performance  of  the  NRS 
filter  with  ramp  failures  induced  on  the  SV’s.  Some  of  the  statistics  are  masked  by  the 
large  vertical  scale  of  the  plots,  but  this  allows  the  full  range  of  the  mean  error  and  residual 
values  to  be  seen,  which  best  characterizes  the  filter  performance.  A  variety  of  GLR  and 
X  plots  illustrate  the  EDI  performance  for  all  ramp  failures  considered.  Several  plots  are 
shown  twice  with  two  different  vertical  scales  to  reveal  significant  details  that  are  lost  in 
other  plots. 
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Figure  G.2,  State  Plots:  Run  8,  Ramp=2T 
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Figure  G.8.  CHI:  Run  8,  Rarnp=2T 


Figure  G.ll.  CIII:  Run  9,  Ramp=lT 


Appendix  H.  GPS  Failure  Plots 


Selected  state  and  residual  plots  are  shown  to  illustrate  the  complete  degradation  of 
the  NRS  filter  with  the  loss  of  GPS  signals. 
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Figure  H.3.  Residual  Plots:  Run  10,  GPS  Fail 
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Appendix  I.  Software 


The  routines  shown  in  this  appendix  were  run  on  Matrix*  (8)  with  the  data  provided 
by  MSOFE  (22)  a,nd  stored  in  Matrix*  compatible  data  files.  The  theory  presented  in 
Chapter  II  is  implemented  with  documentation  to  assist  the  reader, 

//  This  is  a  MatrixX  macro  that  performs  a  OLU  teat  with 
//  n  =  1  for  all  t.  The  routine  "lout"  is  called  to  get  the 
//  the  time-varying  elements  of  F,  H,  gain,  cov,  ft  gamma  stored 
//  in  data  files  font,  hout,  k2-kll,  and  uout. 

//  The  F  matrix  is  discretized  and  windowing  is  applied. 

//  Host  of  the  variable  names  apply  to  the  theory  presented 
//  in  the  thesis  in  Chapter  II  for  the  GUI  test.  Only  10  measurements 
//  arc  included:  SY1-4  and  TR1-6  with  the  baro  anc  velocity  aiding 
//  left  out  as  discussed  in  the  thesis.  NOTE:  This  routine  takeB 
//  between  30-60  minutes  to  run  depending  on  the  machine... be  patient 
// 

glr="exec(’glr’)";  //set  up  macro 
lout.="  exec  (’lout')":  //set  up  macro 

//load  data  -  NOTE:  these  are  commented  out  as  they  loaded  in  already 
//load  ’uout .mxd’ .load  ’hout .mxd’ .load  ’f out .mxd’ , . . . 

//load  ’k2out  .mxd’ .load  ’ki out.  and’ .load  ’k4out.  .mxd’ , . . . 

//load  ’kSout  mrd’.load  ’k6out .mxd’ .load  ’k/out .mxd’ , . . . 

//load  ’k8out .mxd’ .load  ’kUout .Mid’ .load  ’klOout .mxd’ .load  ’kllout.mxd’ ; 

//load  ’simtru.mxd' ; 

// 

//initialize  constant  variables 

// 

118=15;  //ncuher  of  Btates 

siiutime=3599;  //length  of  sinulatiun=3599 

dt»2;  //sample  rate 

ws=1:10;  //set  window  size 

d=r.l;0;0;0;0;0;C;0;0;0]  ;  //failure  matrix 

n=l;  //failure  step  function,  aBsumed=l  for  all  time 

yold*=[0 ; 0 ; 0 ;0 ;0;0; 0 ;0 ;0 ; 0; 0 ;0 ;0 ; 0;0J  ;  //starts  at  0 

b-*0*ones(15) ;  //driving  terms  matrix,  not  really  used=0,  but 

//  needed  for  the  discretize  command 


1-1 


dd=0*ones(10, 15) ;  //D  matrix  not  really  used=0,  but  needed 
//  for  the  discretize  command 
h=0*ones(10,15) ;  //set  up  H-matrix 
f=0*ones(15,15) ;  //set  up  F-matrix 
cov=0*ones(10,10) ;  //set  up  residual  Cov  matrix 
// 

//constant  R  values 

// 

h(l:4,14)=[-l;-i;-l;-l] ; 
h(5: 10,12)=[~1;-1;-1;-1;-1;-1]  ; 

// 

//constant  F  values 

// 

f (9 , 10) =3 . 0668051539128074D-6 ; 
f (10,9)=1; 
f (14,15)=!; 
f (11  .11)=-l/fiOO: 
f (9, 11)=. 0004; 

1(10, 11)=. 03; 

s=0 ; . . .  //s  for  summing 

c.=0 ;  .  .  .  //c  for  summing 

// 

//Initialize  the  window 

// 

lor  k=2:wa-l, . . . 

kt»k,...  //set  index  value  for  lout  routine 
]lout[;...  //get  F,H,gadn(K) ,cov,  l  gamma 
88=[f  b;h  dd] ; . . .  //form  state  apace  matrix 

8fi=discretize(Bs,ns,dt, 'ztransform’)  ;  . . .  //discretize  the  model 
[phi,bd,cd]=split(sd,n3) ; . . .  //split  into  individual  matrices 
covi=inv(cov) ; . . .  //get  residual  cov-inverse 
ynew=phi*(eye(15)-gain*h)*yold-phi*gain*d*n; . . .  //get  y 
yold=ynew; .  . .  //t,tore  recursive  value 
m=h+ynew+d; . . .  //get  m 

sterjn(k)=n’*covi*gamma; . .  .  //find  summation  term  for  s 
8=B-i-steria(k) ; .  . .  //add  sterms  to  get  s 
ctera(k)=si’*covi<>m;  ...  //find  summation  term  for  c 
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c=c+cterm(k) ; , . .  //add  cterms  to  get  c 
l(k)=(s/c)*s; . . .  //MLE  l(k) 
end, 

// 

//Slide  the  window 

// 

for  k=ws : simtime , . . . 

kt“k; .  . .  //get  new  value  in  window 
]  lout  [j .  .  .  //get  F, II, gain (K)  ,cov ,  &  gamma 
s&=[f  b;h  dd]  ; . . .  //form  state  space  matrix 

sd'discretizeCss.ns.dt, ’ztransform') ; . . ,  //discretize  the  model 
[phi  ,bd,cd] “split  (sd, ns)  .  //split  into  individual  matrices 

covi“inv(cov) ;  . . ,  //get  residual  cov-inverse 
ynew=phi*(eye(15)-gain*h)*yold-phi*gain*d*n; . . .  //get  y 
yold“ynew; . . ,  //keep  recursive  value 
m“h*ynew+d; . . .  //get  m 

8terBLCk)=m’*covi*gamiiia; . . ,  //find  summation  term  for  s 
8=8+sterm(k)-sterm(k-W8tl) ; ,  //add  sterms  to  get  s 
ctena(k)**m)*covi*m;  . . .  //find  summation  term  for  c 
c“c+cterm(k)-cterm(k-ws+l) ; . . .  //add  cterms  to  get  c 
l(»  =  (8/c)*0i.  ..  //MLE  1(0 
end, 

// 

//Plot  results 

// 

for  i=l:8imtime, . . . 

time(i)=t(2*i,l)  ;  .  . ,  //keep  time  variable  only  from  simtru.mxd  data 
end, 

plot(time,l, ’repo7/xlabel/Time  (sec)/ylabel/MLE,  L(k)/upper’) 
Wdoorbell  // indicates  end  of  program 

*********************************************  ****************** 
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//Routine  "lout":  Load  in  run  data  for  post  processing.. 

//called  by  macro  "glr". . . 

//... 

//Variable  H,  F,  t  K  values... 

//... 

h(l,l)“hout(kt,l) ; . . .  //H-Matrix 
h(l,2)*>hout(kt,2) ; . . . 
h(l,10)=hout(kt,3) ; . . . 
h(2,l)»hout(kt.4);... 
h(2,2)=hout(kt,5) ; . . . 
h(2,10)*»houtOrt,6)  ; . . . 
h(3,l)=hout(kt,7) ; . . . 
h(3,2)=hout.(kt,8)  ; . . . 
h(3,10)»hout(kt,9) ; . . . 
h(4,l)=hout(kt,10) ; . . . 
h(4,2)=hout (kt ,11) j . . . 
h(4,10)=hout(kt,12) ; , . . 
h(5,l)=hout(kt,.l3)  j .  . . 
h(5,2)=hout(kt, 14) ; . . . 
h(5,10)“hout(kt,15) ; . . . 
h(6 , l)“hout (kt , 16) j . . . 
h(6,2)“hout(kt,17) ; . . . 
h(6,10)=hout(kt,18)  : . . . 
h(7,l)^hout(kt»19)  !••• 
h(7,2)=hout(kt,29) ; . . . 
h(7,l0)=hout(kt,21)j... 
h(8,l)“hout(kt,22)  ; . . . 
h(8,2)=hout(kt,23) ; . . . 
h(8,  lO^houtCkt  ,24)  ; . .  . 
h(9 , l)*hout (kt ,25) ; . . . 
h(9,2)"hout(kt,26) ; . . . 
h(9„10)*=hout(kt,27)  ; . . . 
h(10,l)=hout(kt,28) ; . . . 
h(10,2)“hout(kt,29) ; . . . 
h(10,10)“hout(kt,30) ; .  . . 
f (l,3)“fout(kt,l) ; . . .  //I- Matrix 
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f (l,8)*fout(kt,2) ; 
f (2,3)=fout(kt,3) ; 
f (2,7)“fout(kt,4) ; 
f (3,l)“fout(kt,5) ; 
f(3,2)-fout(kt,6); 

1  (4,2)*:£out(kt,7) ; 
f (4,3)“fout(kt,8) ; 
f(4,S)-fout(kt,9); 

1 (4,6)-fout(kt,10) 

1 (4,8)“fout(kt , 11) 
f (5,l)»fout(kt,12) 

1 (5,3)=fout(kt, 13) 
f  <5 ,4)  «>f  out  (kt ,  14) 
f (5,6)“fout(kt, 15) 
f (5,7)“fout(kt, 16) 
1(6,1) =f out (kt , 17) 

1 (6 , 2) “font (kt , 18) 
f (6 , 4) “f out (kt , 19) 
f (6,5)=fout(kt,20) 
f (7 , l)=fout(kt ,21) 
l(7,2)«fout(Jrt,22) 

1'  (7  ,3>'Tfout(kt  ,23) 
f (T,5)»fout(kt,24) 
f  (7, 6)  “tout  (kt,  ,25] 
f(7,7)«i<wt(kt,26) 
f (7 , 8) out (kt , 27 ) 
f (7, 9) -lout (kt, 28) 
f (8,l)*fout(kt,29) 
f  (8,2)*fout(kt.,30) 
f  (8,3)"fout.(kt,31] 

1  (8,4)-fout(kt.,32) 
f  (8,6)“fout(kt  ,33! 
f  (8,7)  “f  out  (kt ,  34) 
f (8,8)«fout(kt,35) 
f  (8, 9) “lout (kt, 36) 
1(9,1) "lout (kt ,  37 
i(9,2)«fout(kt,38); 


i (9,3)=fout(kt,39) ; . . . 
f (9,4)=fout(kt,40) ; . . . 
f (9,5)=fout(kt,41) ; . . . 
f (9,7)=fout(kc,42) ; . . . 
f (9,8)=fout(kt,43) ; . . . 

gaint(l,l:15)=k2(kt, :)  ; . . .  //Gain  transposed 

gaint (2,1: IS) =k3(kt, 

gaint (3 , 1 :  15)*»k4(kt, :);... 

gaint(4,l : 15)=kb(kt , 

gaint (5, 1 : 15)=k6(kt , 

gaint (6,1: 15) =k7(kt, 

gaint (7,1: 15) =k8(kt, :) ; — 

gaint (8,1: 15) =k9(kt, :);... 

gaint  (9,1: 15)“klOQrt ,:);... 

gaint (10, l:15)«kll(kt, :);.., 

gammatO:  10)=u(kt,  [3  5  7  9  11  13  15  17  19  21]);...  //Residuals 
gamma=gammat ’ ; . . .  //transpose 
gain-gaint ’ ; . . .  // transpose 
ior  j=l:lQ, . . . 

cov(j,j)=sqrt(u(kt,2*j)) . .  //Residual  Cov  diagonal  terms 
end, . . . 


*** ***^4************* t**************^  ****************** ********* 
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//MatrixX  routine  to  generate  a  chi-squared  plot. 

//  Data  files  used  are  uout.mxd  and  simtru.mxd. 

// 

c="exec(’chi,)M;  //  set  up  macro 

clear  r  v  vv  vi  1  n  1  Iw  //  clear  values  from  last  run 

//initialize  values 

k=0; 

nm=10;  //number  of  measurements 

// 

//  Separate  uout  into  r  and  v  vectors.  Where  r  represents 
//  the  residuals,  and  v  represents  the  residual  covariances. 

// 

for  i=3:2:22, . . . 
k=k+l ; . . . 
r( :  ,k)«=u(:  ,i)  ; . . . 
v(: ,k)“u( : , i+1) j . . . 
end 
// 

//  Find  rVr  scalar  for  all  time 

// 

for  j*«2:3599, . . ,  //window  of  time 
vv=0*ones (nm) ; . . .  //initialize  residual  cov  matrix 
for  m=l :nm, . . . 

vv(m,m)=v(j ,m)  ; . . .  //put  diagonal  terms  in  a  square  matrix 
end , . . . 

vi=inv(vv) ; . . .  //Get  vv-inverse 
rvr(j)=r(j  ,  :)*vi#r(j  .  //Gat  summation  term 

time ( j )=t(2* j ,1) j .. .//extract  time  variable  from  simtru.mxd 
end 
// 

//  Find  chi(k)  summation 

// 

chi<l)“rvr(l) ; 
for  k=2 : 3599 , . . . 

chi(k)«chi(k-l)+rvr(k) ; . . . 
end 
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n»15;  //size  of  window 
//  intialize  first  window  values 
for  j“l:n-l, . . . 

lwl(j)=chi(j)  ; . . . 
end 

//  compute  remaining  values  over  time  by  adding  new 
//  values  into  the  window  and  subtracting  old  values  out 
//  of  the  window, 
for  j=n:3599, . . . 

lsl ( j ) =chi ( j ) -chi ( j -n+1 ) ; . . . 
end 

//  clean  up  first  part  of  plot  if  transients  are  too  large 
for  j*=l:300, . . . 

lwi(j)=lwl(3Q0) ; . . . 
end 

n=45;  //size  of  window 
//  intialize  first  window  values 
for  j=l:n-l, . . . 

Iw3(,i)=chi(j) ; . . . 
end 

//  compute  remaining  values  over  time  by  adding  new 
//  values  into  the  window  and  subtracting  old  v  dues  out 
//  of  the  window, 
f or  j =n:3599 j ; t  s 
lw3(j)"chi(j) -chi ( j -n+1 ) ; . . . 
end 

//  clean  up  first  part  of  plot  if  transients  are  too  large 
for  j=l:300, . . . 

Iw3(j)«lw3(300);... 

end 

n“75;  //size  of  window 
//  intialize  first  window  values 
for  , . . . 

Iw5(j)»chi(j)  ; . . . 
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end 

//  compute  remaining  values  over  time  by  adding  new 
//  values  into  the  sindov  and  subtracting  old  values  out 
//of  the  sindos. 
for  j=n:3599, . . . 

Ie5(j)=chi(j)-chi(j-n+l) ; . . . 
end 

//  clean  up  first  part  of  plot  if  transients  axe  too  large 
for  j=l : 300 .... 

Iw5(j)“ls5(300) ; . . . 
end 

Wcuckoo  //  indicates  end  of  routine 

//plot  results 

plotftime, [1b1,1«3,1w5] , ’repo7/xlabel/Time  (sec)/ylabel/Chi(k)/. . 

yiuiu-'G/iegeuu/io  V juiuO w  |  iG  VliiidOi  !  7G  window/  . 
title//') 
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